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ON THE SUBNORMAL VARIANCE IN THE COUNTS OF RAN- 
DOMLY DISTRIBUTED PARTICLES: I. APPROXIMATE 
TREATMENT OF THE THREE-DIMENSIONAL CASE AND 
DISCUSSION OF EXPERIMENTS pra 


HESSEL DE VRIES* 


RiyKs UNIVERSITY 
GRONINGEN, NETHERLANDS 


The question is raised concerning the possible causes of abnormally small standard devia- 
tions found in counting samples in which particles are distributed at random (e.g., blood cells, 
fat globules in milk, etc.). The effect of discarding abnormal samples is discounted inasmuch as 
small standard deviations occur even when all samples are counted. An approximation method 
is used to calculate the effect of finite particle size, of known repulsive forces between particles 
and of convection currents. This calculation shows that neither finite size nor the known 
repulsive forces are sufficient to account for the observed abnormality of standard deviation, 
but that convection currents can possibly account for it. The possible presence of long-range 
repulsive forces cannot, however, be excluded. 


1. Introduction. As early as 1922, R. A. Fisher, H. G. Thornton, and 
W. A. Mackenzie (Fisher, 1950) discussed a case of subnormal variance in 
the case of bacteria counts. J. Bergson et al. (1940) have reported a vari- 
ance of 0.92 times the normal in blood cell counts, and A. van Kreveld 
(1938, 1942, 1946, 1947) has discussed careful counts of bacteria and of fat 
globules in milk, where again a subnormal variance was found. The pres- 
ent study started from another case of subnormal variance in blood cell 
counts, reported to Dr. A. Rapoport by Dr. Gerda Seidelin of Copen- 
hagen. 

No final explanation has been arrived at concerning these persistent 
occurrences of subnormal variance. The deviations discussed by Fisher et 
al. were, at least in one case, related in some way to an error in the experi- 
mental procedure, and from this it was concluded that “subnormal vari- 
ance is a danger signal which cannot be disregarded.” Indeed it is easy to 
see how a wrong procedure can lead to subnormal variations in making 
blood counts. An experienced observer will immediately see that some of 
the squares contain an “abnormal” number of cells and will not waste 

* The main part of this work was done in the Committee on Mathematical Biology at the 
University of Chicago and originated with a suggestion by Dr. A. Rapoport. 
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time counting them. Moreover if in counting a sample he sees that the 
result deviates too much from the approximate mean, he may discard the 
sample. But such explanations would not apply to the counts reported by 
Bergson et al. and by van Kreveld since those counts were made with the 
specific purpose to determine the standard deviations in randomly counted 
samples. Bergson reports the suggestion that the finite dimensions of the 
particles counted might account for the subnormal variance. To a certain 
extent this also covers the suggestion of Dr. Seidelin* that repulsive forces 
might play a part. Apparently both of these factors would tend to dimin- 
ish the possibilities for fluctuation. In fact, strong repulsive forces or large 
dimensions resulting in nearly complete filling of the available space 
would give rise to a very regular distribution of the particles. In the case 
of bacteria, repulsive forces could be stimulated by chemical agents which 
could have an inhibitive effect on growth or a repulsive effect. These 
explanations will be discussed in more detail in the following sections 
together with the thorough study by van Kreveld (loc. cit.). 

A. Rapoport pointed out that the problem is analogous to that of a 
neuron responding to a Poisson-distributed shower of stimuli, where the 
regularity of the response increases with the length of the refractory pe- 
riod. This problem also arises in nuclear physics where the finite duration 
of a response of a counter corresponds to the finite geometric dimensions 
of the particles in our problem. C. Levert and A. L. Scheen (1943) and 
L. Kosten (1943) have solved the counter problem. However, they con- 
sidered a counter which was ready to respond only when, during its refrac- 
tory period 7, no new particles came in. Our problem is more similar to a 
counter which is ready to respond to the first particle which comes after a 
time 7 after the last particle which was recorded. This case was considered 
by H. A. Van de Velden and P. M. Endt (1942). However, they gave only 
an approximate solution. Finally one meets the same problem in the the- 
ory of the distribution of molecules in gases. 

2. Elementary calculation of the effect of finite dimensions. There are two 
ways of calculating the standard deviation o of the particle counts. The 
first method starts with deriving the distribution law which replaces the 
well-known Poisson law if the assumption of finite size is made. We suc- 
ceded in deriving this law, by utilizing some expressions derived by F. Zer- 
nike and J. A. Prins (1927). This new law proved so difficult to handle, 
however, that subsequent calculation of ¢ was not undertaken. 

We shall therefore proceed otherwise. Suppose vanishingly small linear 
particles are distributed over a line segment of length L. We suppose this 


* Private communication to Dr. Rapoport. 
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line segment divided into N equal intervals, V being so large that one 
interval will not contain more than one particle. Suppose we repeat our 
process of randomly distributing the particles over similar line segments, 
keeping the over-all density constant so that the mean number of particles 
on L is x. We denote by x; = 1 the presence of a particle in the &th in- 


terval and by x; = 0 an absence. Then for any interval we have 


oa a— (z)*= > (1-4). (1) 


Let y,, denote the total number of particles in the first m intervals. Adding 
one interval to these m intervals, we have: 


Oa oe te, te Yad ae (2) 
If the particles are vanishingly small, there is no correlation between 
%m4i and the situation in the other intervals. In that case, the expression 
in parenthesis of (2) vanishes, and we obtain by induction the well-known 
formula of Bernouilli 


oi, = Not=ni( 1-4). (3) 


As WN becomes definitely large, this reduces to 0? = n. 

Suppose now the particles have a finite length a. Then there is a cor- 
relation between ¥mii and y» in (2), so that m4 no longer equals 7 nx. 
The product yn%m41 Vanishes unless #41 = 1. Then, however, the mean 
value of y is smaller than mn/N, since the last intervals to the left of #41 
have to be empty (we mark the position of each particle by its left edge). 
Our problem is now to determine the new value of jm, which we will de- 
note by Fn. 

Let f = an/L stand for the fraction of the line which is occupied by the 
particles. Then if f is small, we can say, as a first approximation, that the 
length of the line to the left of x,,.41 is reduced by a, that is, the length of 
one particle. Then we have for the value of 5m the expression (m 
— aN /L)n/L, instead of mn/N. Substituting this value into equation (2), 
we find for the whole line and large NV 


ot=n(l—2f).. 4) 


A similar calculation can be given for the two- and three-dimensional 
cases if the particles are considered circular or rectangular. The total space 
can be divided into small parts appropriately numbered. When the 
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(m + 1)th part is occupied by a particle, one will find that of the m fore- 
going parts an area fwice the area of one particle has to be empty. In the 
three-dimensional case, the empty volume turns out to be four times the 
volume of one particle. Thus the correction factors corresponding to 
(1 — 2f) in equation (4) turn out to be (1 — 4f) and (1 — 8f) for the 
two- and three-dimensional cases respectively. The formulae hold ob- 
viously only for small values of f, since for f = 1 (space completely filled 
with particles) the fluctuations should vanish, whereas the formulae yield 
negative values for sufficiently large f. 

The error arises from the assumption that 9, corresponds to the mean 
number of particles on a line which is reduced in length by a. This shorter 
line, however, is not a line randomly chosen but adjacent to its right-hand 
end there is an empty space, so that we estimate 7 too small. For small 
values of f, the difference will be negligible, since a line chosen at random 
would also have empty adjacent regions in most cases. 

For the one-dimensional case, we can compute 7,, exactly. If the line is 
chosen at random, its right end will either fall on a particle (with probabil- 
ity f) or on an empty space. In the second case, it will contain 7’ particles 
according to our definition of 7’. In the first case, the particle will project 
to the left over a distance, whose expected value is a/2, so that the mean 
number of particles will be 7’(Z — a/2) + 1. Now y'(Z — a/2) will equal 
y’L — f/2, since f/2 is the mean number of particles ona length ia. 
Strictly speaking, this is true only when L is so long that there is no cor- 
relation between the positions of the particles at the two ends of the line 
(Zernike and Prins, 1927). 

The smaller the empty spaces, the larger should L be for this inde- 
pendence between the two ends to hold.* With this restriction, we find 


gy=s'+s(1-4), (5) 


which expresses j as the weighted sum of the values of y’ in the two cases 
considered. Substituting this result into (2), one obtains 


ae Ca) ae (6) 


which agrees with (4) for f <1. For larger values of f, formula (6) should 
be used. 


* For f = 1, there will be a correlation between the ends of the line whatever the value of 
L. For f = 1, however, the calculation is very simple. If we write L = pa + Ba, where p is an 
integer and 0 < 6 <1, then o* turns out to be (1 — 8). Unless % is small or f very nearly 
unity, this can be neglected relative to the fluctuations given by our formula. 
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An analogous exact treatment of the two- and three-dimensional cases 
is much more complicated. The three-dimensional case has been consid- 
ered in detail by van Kreveld (1938). He also arrives at the formula 


at 2 ais (7) 


as a special case of more difficult expressions, which were not pursued 
further. Since only this special case could be computed, it seemed worth 
while to give the above elementary approach. van Kreveld’s reference to 
the kinetic theory of gases, however, opened the possibility of calculating 
higher order terms, which became important when f is no longer very 
small compared to unity. 

3. The influence of repulsive forces. The only repulsive forces between 
suspended particles which are known from the theory of colloids are elec- 
trical forces, or, rather, the forces between the ionic double layers around 
the particles. Their influence would be to give the particles a larger “‘effec- 
tive volume.” As an approximation, we can say that the new boundaries of 
the cells should be chosen so that the potential energy equals kT when the 
enlarged cells just touch one another,* where & is Boltzmann’s constant 
and T the absolute temperature. Formulae for the potential energy per 
cm.” of two parallel planes have been given by E. J. W. Verwey and J. Th. S. 
Overbeek (1948). Different fluids are used for diluting the blood in the 
counting chambers, but all of them would give less repulsion than a 1% 
NaCl solution (the forces decrease with increasing concentration). Start- 
ing from these considerations and the formulae given by Verwey and 
Overbeek (cf. Table XI, pp. 83 and 85, loc. cit.), we find that the thickness 
of the ionic layer which has to be taken into account is about 10-° cm. 
This is the estimate if the blood cells are considered to be in parallel planes 
just opposite each other. For other relative positions, the correction is 
even smaller. Thus the correction is of negligible proportions, as might be 
expected from what is known about these double layers from other crude 
calculations. It appears, therefore, that electrical repulsive forces cannot 
play an important part in influencing the distribution of particles. 

4. The effect of finite volume in blood counts. The average diameter of 
blood cells is about 9 w, the thickness is 1.5 u. The counting chamber 
measures 200 X 200 X 100 uv? and contains on the average about 100 


* van Kreveld (Joc. cit.) also gives formulae for the influence of repulsive forces on o?. The 
effective volume which emerges from his reasoning is practically the same which we assume. 
Since it is shown here that the effect of the electrical forces is not important, we will not go into 


more detail. 
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cells. So the fraction f of the volume occupied by the cells is 23 X 1052: 
We find 
o=n(1—8f) =n(1—0.018) (8) 


o=0.99Vn. 


Hence the correction factor amounts to 0.99 and cannot account for any 
appreciable deviations from the normal law. 

It has to be pointed out, however, that the blood cells are generally 
counted after they have settled on the bottom of the counting chamber. 
So the final density and hence f is much larger than what we calculated 
above. Of course, this effect cannot change the total number of cells in the 
chamber as a whole, but in fact what was called a counting chamber is 
only a little square in the larger hole in the object glass into which the 
droplet of blood is brought. This remark could also be important for the 
counting of fat globules in milk. Here the fat globules move upward to the 
cover glass, and again only a small fraction of the total area is counted. 
Consequently, also here a regular distribution of the particles over the 
whole cell would reduce the fluctuations per counted area below the nor- 
mal value. In part of the counts of fat globules, however, the particles 
were “‘fixed”’ by gelatine, so that they could not move upward and gather 
on the cover glass. Even here, however (van Kreveld, 1942), a subnormal 
variance was found. 

5. The influence of convection currents. Suppose that in a certain area 
the number of blood cells is larger than normal. Then the specific weight 
will be larger here, so that convection currents will flow which tend to 
remove the variations. New irregularities, however, can be brought about 
by diffusion. Thus we see that a counter action exists between convection 
currents and diffusion and that the “normal” fluctuations could be partly 
suppressed by convection currents. 

An exact calculation of this effect is very difficult, but an estimate can 
be made as follows. We consider two cubes, each of side x, which contain 
different numbers of particles and consequently have different weights. 
Then the cubes will move relative to each other. We wish to calculate the 
time ¢., after which the displacement is x. Though this does not yet mean 
complete mixing, it is a good indication of the velocity of the convection 
current. (The current need not be exactly vertical, but the heavier part 
may flow partly under the lighter part.) One finds ¢, = «?n/K, where K is 
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the force and 7 the viscosity of the fluid which we will suppose to be water 
(» = 0.01). The force K is given by 


K=.05«3 gf AC dynes. (9) 


Here g is the acceleration of gravity, f the fraction of the fluid occupied by 
particles, and .05 is approximately the difference of the specific weights of 
fluid and particles (for blood cells as well as fat globules in milk). In both 
cases, f is about 0.2% (cf. Section 4 above). Finally AC is the difference in 
concentration of the particles in the two cubes, relative to the mean con- 
centration C. Approximately, (AC)? will be inversely proportional to the 
number of particles, and therefore to x°. The factor of proportionality 
follows from considerations given in Section 4 (in a volume of 4 X 10-3 
mm.*, AC was 0.1). For fat globules in the same volume of milk (also 
diluted before counting), AC is somewhat smaller (cf. van Kreveld, loc. 
cit.). We finally find 


i.v500 Vx. (10) 


It follows from (10) that convection takes more time in larger volumes. 
Convection is of interest for our problem only when it ranges over at least 
the dimensions of the counting chamber, i.e., about 0.2 cm. (cf. above). 
Substituting this value for x in (10), we get ¢, = 70 sec. Since the whole 
procedure of counting, etc., takes more than 70 sec., we conclude that 
convection could have some influence on the distribution of particles so as 
to reduce the normal fluctuations. 

This calculation, however, must be supplemented by an estimate of the 
counter action of diffusion, which “‘tries’”’ to restore the normal random 
distribution. From the well-known formula of Einstein for Brownian mo- 
tion of particles, we find 

ta = 10°x?. (11) 


Here x is the displacement of the particles which we suppose to be 
spheres with radius 6 y. Inserting here also « = .02, we find it would take 
about 4 X 10° sec. to reach the normal random distribution. Comparing 
(10) with (11) we see that even in the small volume of a counting chamber 
convection can be considerably faster than diffusion. It is only a crude 
estimate, however, which at most indicates that further work in this direc- 
tion may be useful. 

Even the convection effect, however, may not explain all subnormal 
variances reported. In a discussion with Dr. van Kreveld, he drew the 
author’s attention to the marked influence of various chemicals on the 
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fluctuations in the counts of fat globules (van Kreveld, 1946). It is difficult 
to see how these chemicals could have an appreciable influence on convec- 
tion. Therefore van Kreveld still supposes that long range repulsive forces 
may exist, though their nature is still unknown. Further experimental 
work seems to be indicated. 
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RANDOMLY DISTRIBUTED PARTICLES: II. EXACT 
TREATMENT OF THE ONE-DIMENSIONAL CASE 
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The effect of finite particle size on the standard deviation in sample counts is computed for 
the one-dimensional case. To a first order of approximation the correction is found to be identi- 
cal with that found by H. de Vries (1953) using a general approximation method. 


H. de Vries (Joc. cit.) used an approximation method to derive the cor- 
rection factors which must be applied to the variance in the density of 
particles randomly distributed in one-, two-, and three-dimensional space 
when the fraction of the space occupied by the particles is finite. We shall 
consider this problem further from the theoretical point of view. In par- 
ticular we shall rigorously derive the result for the one-dimensional case. 

Consider a straight line (which may, for example, represent the time 
axis) upon which an infinite number of one-dimensional rods (or “par- 
ticles”) of length a are distributed at random but in such a way that no 
two rods overlap. If / is the average separation between the particles, then 


(1) 


a 


I+a 


represents the fraction of the line occupied by the particles. The number 
of particles within a certain interval of length «—or, more accurately, the 
number of left-end points within that interval—will be denoted by n(x). 
Its average value, when the interval of length x is chosen in an arbitrary 
way, is evidently given by 


f= 


x 
Zz 
ak (2) 


We are interested in the fluctuations of (x), that is, we shall compute 
the quantity 


22) 


o?(x) =n? (x) — [n(a)]? nis (3) 


for an arbitrary length x. 
7 245 
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In what follows we shall take / + a as unit of length, so that the number 
density, that is, the average number of particles per unit length, is taken 
to be one. We start with the radial distribution function g(x) in one dimen- 
sion, which was first introduced by F. Zernike and J. A. Prins (1927), 
defined as follows: g(x)dx measures the probability of finding the left-end 
point of a particle in the interval (x, « + dx) for positive « when the left- 
end point of one particle is known to be situated at x = 0. Zernike and 
Prins have shown that 


g(x) =0 forx<a, 
x—ka 


tee 7 pqvtci P| (4) 
pe) ae ra>e.| 


The kth term in the sum gives the probability of finding the &th particle 
in the interval (x, « + dx). It will be seen that g(x) is not an analytic 
function, even for x > a, since the (7 — 1)th derivative of g(x) is dis- 
continuous at the points x = na, where x is an integer. It turns out, how- 
ever, that the Laplace transform of g(x), defined by 


Z(s) =f ep (x) da 
0 


has a very simple form. In fact it can be easily shown that 


s be = e—kse Sy 1 
ats) 2 GED, Gl+i) eat 6) 
Since for large s, g(s) = (sle**)—", we see that 
lim g(x) = Je 
xa 1 
as can also be seen directly from (4). At s = 0, 2(s) has a singularity like 
that of 1/s. From this it follows with the help of one of the Tauber 


theorems (Doetsch, 1937) that for large « g(x) approaches unity ex- 
ponentially. Furthermore, writing 


we may expand g(s) and therefore also g(x) into powers of f, as is also cus- 
tomary in the three-dimensional case where expansion is made in terms of 
the powers of the density (de Boer, 1948—49; Nijboer and Van Hove, 
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1952). Thus 


g(x) = >> gn (x) fm. (6) 
n=0 


It can be shown that for 2 > 0, g,(x) = 0 for x > (n+ 1)a. Finally 
we notice that for / = 0, 


Hence 


g (x) = S96 (x—na) , 
eI 


where 6(«) represents Dirac’s 6-function, so that we have a rigid one-di- 
mensional lattice, as, of course, should be the case if the line is completely 
filled with the particles. 

Let us now call f,(«) the probability of finding m (left-end points of) 
particles in an arbitrarily chosen interval of length x. Obviously there is 
a close connection between the functions f,(~) and g(x). In particular, 
their Laplace transforms are simply related. One can show that 


Jn (s) = SPA SAE (st 1) om)? (7) 


If we put a = 0, we find 
ar — B 
lim fn (Ss) ~ CsI 1)? 
so that 
exp (—4) 
2 ty LA l 


which is Poisson’s formula, as should be the case. Also in the limit of 
vanishing /, we find a simple result for f,(«). 

The average values of » and n? can now be easily calculated in the 
image space of the Laplace transform. In fact, we have 


ae 1 
en tals) are (8) 


n=1 


so that n(x) = x, as we expected, and 


a RRR oO cciabi ie 
Dyn (m2) Inls) SSE oi (9) 
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Comparison of (5) and (9) shows that o? may be expressed as a convolu- 
tion, namely, 
ot=2 f (x—1) gO dita—a?. (10) 
0 


Formula (10) is particularly useful for calculating the fluctuation for small 
values of x, because in this case g(x) consists of a small number of terms 
only, as can be seen from (4). In particular for « < a, we have g(x) = 0, 
so that 
o?=x—x? for <5 

In order to find o? for large values of x, one may, according to one of the 
well-known Tauber theorems (Doetsch, 1937), investigate the behavior 
of its Laplace transform in the neighborhood of its singularity s = 0. 
From (8) and (9) we derive the following expression, 


: lee" 
a PH 2aes+e 
G(s) soe AE (11) 
-+ power series in s. 
Therefore / : 
e(x) =a(1— fet 1—-Z4s (12) 


+ remainder term. 

The remainder term vanishes exponentially for large x [except for f = 1, 
but in this case the fluctuations can be found directly from (9)]. From 
(10) it follows easily that the remainder term has the form 


—2 f° (e- {g( -1) at. (13) 


If f «1, the remainder term can with the help of (6) be expanded 
into powers of f. The first term of this expansion is the term with f* if 
x > (n+ I)a. 


For infinitely large « and f # 1, we see from (12) that 


gt mee fo LE Le (14) 


which is the result obtained by de Vries (Joc. cit.) by elementary methods. 
Equation (14) holds approximately for finite large x, the better the 
smaller f is. Indeed, closer examination shows that (14) is a good approxi- 
mation if only « > a(1 — f)~. 

If one is interested only in the fluctuations in the number of particles in 
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an infinitely large interval (or volume), then the easiest way of deriving 
formula (14) is from the equation of state of the fluid formed by the par- 
ticles considered. According to a well-known result due to Gibbs, one has 
(de Boer, 1948-49) 


FS hy (3 


n (co) ge (15) 


The left side represents the relative fluctuation in a large test volume 
(which, of course, should still be compared to the total volume). The right 
side, where v denotes the volume per particle, p the pressure, T the abso- 
lute temperature, and k Boltzmann’s constant represents the relative com- 
pressibility of the gas considered, that is, the compressibility of the gas 
divided by the compressibility of an ideal gas of the same density. By 
means of this formula one could also take into account attractive or repul- 
sive forces between the particles, if they exist. 

Let us first apply (15) to the case considered above, which translated 
into kinetic theory terms denotes a one-dimensional gas of hard spheres. 
In this case, the equation of state is simply (de Boer, 1948-49) 


p(v—a) =hkT, 
so that 


ats PaaS a 
n (o) vp? 


“y= a-/" (16) 
as is given by (14). 

The determination of the radial distribution function for a three- 
dimensional gas of hard spheres seems to present considerable difficulties 
(Nijboer and Van Hove, 1952). We shall consider only the fluctuations in 
large test volumes. Then we can apply (15) again. For small densities, the 
equation of state can be expanded into powers of the density. Only the 
first few coefficients in this expansion are known for a gas of hard spheres 
(Fowler and Guggenheim, 1939; Nijboer and Van Hove, 1952). They are 
given by 


2 53 
pom a(1454+3 540.2869 54...), (17) 


where » is the volume per particle and 0/4 is the volume of one particle. 
For higher densities, approximate expressions for the equation of state 
have been given, but will not be discussed here. ; 
Substituting (17) into (15), we have 


eee Sy toidy? = LO Ge itn, (18) 
n (@) 
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where f = b/4y is, as before, the fraction of the total volume occupied by 
the particles. If f is very small, one can be satisfied by taking the first two 
terms only, which is again a result given by de Vries (Joc. cit.). 

In connection with a problem closely related to the one treated here, 
van Kreveld (1947) had already applied formula (15) to the case of a van 
der Waals gas. 
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A critical examination of the ‘“‘classical” theories of photoreception in view of more recent 
experimental findings yields the result that these theories do not possess the property to de- 
scribe all the more significant phenomena of photoreception correctly, and to some extent suffer 
the lack of more general applicability. The basis for a new and presumably more general theory 
of photoreception based on dynamical aspects is laid out. Emphasis is put on the time course 
of afferent and efferent excitation in the photoreception model, consisting of a receptor element, 
an afferent and an efferent neuron of the one-factor Rashevsky-type, and an effector organ. 


I. Introduction. The process of photoreception with its major phenome- 
na such as light susception, latent and reaction times, light and dark 
adaptation, intensity discrimination and visual acuity has long been an 
object not only of numerous extensive experimental studies but also of 
some theoretical ones.* Among the latter the most thorough and extensive 
one, that of S. Hecht (1931, 1934), will be briefly considered here. 

In his photoreception model consisting of a receptor element only, 
Hecht assumes two groups of substances located in the receptor. One 
group consists of m photosensitive substances S = 5, ... Sm, which react 
under the action of light to produce a second group P = fi... pn of n 
reaction products. The velocity of the “light process,” S — P, is propor- 
tional to the intensity of light, J, and the concentration of S, [S] = 
(a —- x). Simultaneously and independent of light the restitution of the 
photosensitive substances, P — S, takes place in a reversible process with 
a velocity proportional to the concentration of P, [P] = x. Without any 
loss of generality we may assume a = 1, so that the whole process, 


SP, follows the differential equation 


CART (feeeey aes bia (1) 
dt 
* The conception of photoreception as adopted here will be confined to the phenomena men- 
tioned, i.e., to the more primary process of vision not necessitating a complicated poe 
much more complicated, secondary processes such as perception of space and space-time differ- 
ences (gestalt, movement, etc.) should be referred to as photoperception. 
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which integrates for m= 1 and n= 2 (claimed adequate for visual 
purple by Hecht) to the form 
1 -- Cy | 
mr pete el ee, Se (2) 
e Ap ptearer : 

with a, c, and c, as constants for x(t = 0) = 0 and J = constant. Now 
according to Hecht’s theory a reaction occurs when x exceeds a threshold 
value. This, of course, implies that the excitation intensity in the receptor 
element is a function of x. Hecht obtains, for J = constant, an excitation 
intensity in the receptor element as a monotonically increasing asymptotic 
function of time (Fig. 1, curve 1). In view of the experimental findings of 


t 


FicureE 1. Time course of excitation intensity. Curve J: after Hecht, equation (2); curve 2: 
after Hecht, with a one-factor neuron; curve 3: idealized action potential P/J after Granit and 
Riddell (1934). For the reader’s convenience the slopes, practically linear in the range of small 
t, of the three curves have been chosen different. Keeping this in mind, it can be seen that 
curves 1 and 2 approximately correspond with curve 3 for the range of small and very large /, 
but not for the intermediate range. 


R. Granit and H. A. Riddell (1934) and others, this is incorrect since the 
action potential, i.e., the excitation intensity of the afferent neuron, de- 
scribes a course having a maximum (Fig. 1, curve 3). This inadequacy of 
Hecht’s theory is not to be ascribed to the lack of an afferent neuron, 
since the introduction of a one-factor type neuron stimulated by a quan- 
tity proportional to « will not yield basically different results (Fig. 1, 
curve 2), As might be seen from a more intense study of Hecht’s works, his 
otherwise excellent theory applies only to the limited range considered by 
him: the very start of the excitation course and the steady state in the 
asymptotic region. Though the application of Rashevsky’s two-factor the- 
ory of nerve excitation and inhibition (Rashevsky, 1948, chap. xxiv) 
might yield satisfactory results (working on Hecht’s basic assumptions), 
another theory based on to some extent new, presumably quite plausible 
dynamical aspects will be developed presently. 

Brief attention has to be drawn to R. Runge’s theory of photosensitiv- 
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ity of lower animals (1945). He assumes direct stimulation of nerve fibers 
by light. The correctness of his assumption will not be discussed here; but, 
on the other hand, his using a direct influence of the stimulus light on the 
threshold value seems to be quite questionable. In addition to statistical 
variations, fluctuations in threshold values under the influence of a stimu- 
lus are established experimentally, but those most probably will have to be 
ascribed to causes other than a direct effect of light on the threshold. Fur- 
thermore the sigmoid time course of excitation, corresponding to Runge’s 


mfe)——-—~< 


Ficure 2. Photoreception model, consisting of receptor element RE, afferent neuron N b 
efferent neuron Ns, and effector organ M. 


theory, does not agree with the experimental results. Runge, at any rate, 
does not claim generality for his theory; thus it is not worth while to 
argue these points. 

On the whole it must be acknowledged that there seems to be no theory 
of photoreception in existence which describes the major phenomena of 
photoreception correctly along with the experimentally established time 
course of neural excitation. 

IT. The photoreception model. The model adopted for the purpose of 
these investigations essentially consists of three parts (Fig. 2): A receptor 
element RE which, upon application of the stimulus light energy, trans- 
forms the absorbed energy into receptor energy. The receptor energy acts 
upon an afferent neuron N, as stimulus with the result that 1; itself be- 
comes excited. This excitation is transmitted to a synapse where, by action 
of the afferent excitation as stimulus, an efferent neuron N» will be excited 
as soon as the afferent excitation intensity exceeds a constant threshold 
value. At the end of NV, an effector organ M will react on the efferent ex- 
citation as soon as the intensity of the latter exceeds a constant threshold. 

1. The receptor-element RE. In the receptor element the primary pho- 
toreception process, the absorption of light energy and its transformation 
to receptor energy, takes place. A light-absorbing substance, or group of 
substances, is assumed to be present. Now, as for the behavior of this sub- 
stance under the action of light, two alternative possibilities seem to exist: 

Case (a): The substance, under the action of light, remains unaltered; 
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it serves as a means of transmitting the light energy only, the latter acting 
thus on a photosensitive substance as in Case (b) below, or on nervous 
fibers directly. We will call this case photochemical sensitization. We shall 
not discuss here whether such a system might be realized in the animal 
kingdom, particularly since it has no bearing in principle on the presumed 
generality of the formulations, as will be shown presently. 

Case (b): Under the action of light the light-absorbing photosensitive 
substance, or group of substances, S, will be converted into a reaction 
product, or products, P, with a velocity proportional to the absorbed light 
energy and the concentration of S. We shall call this process photochemical 
reaction. 

Since, for the sake of simplicity, no passing of substances S or P to or 
from the receptor element will be assumed, the restitution of S will be 
furnished by a reversal process, P — S, with a velocity proportional to the 
concentration of P, independent of light. As for Case (a) of photochemical 
sensitization, the photosensitizing substance can be neglected altogether 
if its task is the transmission of light energy to a second, photosensitive 
substance only. On the other hand, for the direct transmission of light 
energy to nervous elements by means of such a substance one only has to 
hold the concentrations of S and P constant in order to achieve the cor- 
responding results. Therefore Case (a) is a limiting case of Case (6), and 
we may limit ourselves to an investigation of (0). 

So far Hecht’s model has been adopted. But his view that the concen- 
tration of P represents the excitation, or rather acts as a stimulus, has to 
be abandoned. It seems implausible that a process as continuous as nerv- 
ous excitation should be kept flowing by the mere presence of a substance 
in a low, finite concentration. The failure of Hecht’s theory to represent 
the time course of excitation adequately may well be attributed to this 
shortcoming. 

The photochemical reaction, which we suppose for simplicity to be the 
conversion of one photosensitive substance S to its reaction product P 
under the action of light, and the re-synthesis from P 

light 
S.——— P (3) 
“dark” 
may have either one of the following thermodynamical properties: 

(1) The process, S — P, is endothermal, i.e., S is exothermic. The light 
supplies the energy essential for the reaction of the exothermic S to the 
endothermic P; while the reverse reaction P — S is exothermal and spon- 
taneous. In this case (3) would be a reversible process. 


. 
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(2) The process S — P is exothermal, i.e., now S is endothermic. Conse- 
quently, this process would be spontaneous would it not be for its pre- 
sumably having to pass through an energy level higher than that of either 
S or P. Thus this process represents a “lock mechanism,” a mechanism 
quite plausible for biological systems. By absorbing light energy, 5 is ac- 
tivated to react into P, releasing at the same time a possibly much greater 
amount of energy than the activating light energy supplies. A chain reac- 
tion, however, is prevented, as will appear below. The energy for the re- 
verse process, P — S, which is endothermal, has to be supplied from the 
“outside.” Assuming for the source of this energy an enzymatic process, 
€.g., a “redox’’ system, is not too implausible. The process (3), on the 
whole, would then be a pseudo-reversible one. 

Keeping in mind that the object of the receptor element is the transfor- 
mation of the absorbed light energy into receptor energy, the question is: 
What is the nature of this transformation? Taking it for granted that this 
energy will be of a chemical or electrical nature which, under biological 
aspects, would simplify matters enormously, one has in case (1) to assume 
that receptor energy be proportional to the energy set free by the exo- 
thermal process P — S. This assumption was used at first but had to be 
abandoned because it did not yield satisfactory results. On the other 
hand, considering that quite a bit more light quanta per unit time will be 
absorbed than will be needed for the conversion of S to P, one may safely 
assume that the surplus light quanta absorbed, or a proportional quantity, 
represent the receptor energy. However in this case no explanation is made 
of how this quanta shall act as a stimulus on a nerve fiber. Whatever its 
limitations we will adopt the following assumption based on the model 
given in (2) above. 

The receptor energy will be a proportional quantity of the energy released in 
the exothermal process S — P. Transmitted as a stimulus to a nerve fiber this 
receptor energy is “consumed,” thus insuring that no chain reaction will occur 
in the exothermal process. 

2. The afferent neuron N;. No threshold value will be assumed for the 
stimulus receptor energy on the afferent neuron, because the receptor 
element might well be considered as a part of the neuron itself. Of course 
this would imply that NV; is the neuron of the anatomist rather than the 
neuron, or pathway, of the mathematical biophysicist; but neither view 
will be of principal consequence. Owing to the stimulus receptor energy 
the neuron will become excited, the excitation resulting in the formation 
of an excitatory factor at the synapse with the efferent neuron. 

3. The efferent neuron >. This will become excited as soon as the excita- 
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tory factor of Ni equals or exceeds a threshold value, eliciting a reaction 
of the effector organ M as soon as the excitatory factor at the connection 
of N. with M equals or exceeds a threshold value characteristic of M. 

Both neurons, NV, and N2, shall be of the one-factor Rashevsky-type. 

III. The process of photoreception. Summing up, the postulates of the 
theory are: 

1) In a receptor-element RE two substances are present: An endo- 
thermic substance S, which under the action of the activating light energy 
is converted into the exothermic substance P with a velocity % propor- 
tional to the absorbed light energy J’ (the intensity of light falling on the 
receptor element being J) and the concentration of S, [.S]. The necessary 
supply of S is insured by a reversal process, where S is resynthetized 
enzymatically from P with a velocity v2 proportional to the concentration 
of P, [P] = c; no passing of either substances to or from RE is to occur. 
The velocity of the over-all process (3) is then V = 2; — w. 

2) Owing to the exothermal character of the process S — P, an amount 
of energy per each “molecule” of S reacting into P is set free under the 
activating energy of light. This energy, or rather a proportional amount of 
it, is the receptor energy. 

3) An afferent neuron J, is connected with the RE. The receptor en- 
ergy serves as a stimulus to NV; which becomes excited, resulting in the 
production of an excitatory factor e at the synapse with an efferent neuron 
N2 with the threshold . The neuron NV, becomes excited when the excita- 
tory factor ¢ is equal or greater than h. As soon as the excitatory factor ¢’ 
of 2 exceeds the constant threshold h’ of an effector organ M, a reaction is 
elicited (Fig. 4). 

Taking J as the intensity of light falling on RE, k’ the absorption con- 
stant of S, k a reaction constant, and the concentration of S as [S] = 
1 — [P] = 1 —<c, the velocity of the process S — P will be 


11 = kJ+ [1— e709] (4) 
or, for small values of k’(1 — c) as a first approximation, 
1 = kJ (1— ¢) (5) 


with kk’ = ky. Correspondingly, the velocity of the reverse process, 
P — S, will be 
vo = kec. (6) 


Hence the velocity of the over-all process (3) is 


d 
Vas = kJ (1— 0c) kee. (7) 
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Integrating (7) for c and substituting in (5) with c(t =,0) = cs) and 

J = constant 
Vy eee Mit 52 {Ro+ [RiJ — co (ki + ko) ] e ~@sI+k,)t 

hel + bs 2 1 0 1 2 € ; : Be (8) 

Now the receptor energy (playing the role of a stimulus) released from 

RE per unit time element is, of course, proportional to 2, and the inten- 

sity of the stimulus acting on Nj is o = an, a monotonically decreasing 

function of time. That means that just at the very first stage in the photo- 

reception process a “‘primary adaptation” takes place. Though for very 
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meee atlas eT ees 
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Ficure 3. Primary adaptation curves. Time course of the conversion of S to P under the 
action of light in the dark adapted system, for two different light intensities. Equation (8). 


large values of J, o tends to infinity for small values of ¢ and for co < 1, 
it assumes an asymptotic value for t + ~. The corresponding asymptotic 


value of v,; then amounts to 
Ry Rot 
=_ 9 
eS EN, se 


which for J > © attains a finite maximum 2%, — ’». Figure 3 shows the 
“primary adaptation curves” 2 = f(t), where the circumflex (a) eas 
throughout this work, denotes that the initial conditions are taken zero, 
for two different light intensities. For the production of the excitatory 
factor € at the synapse V; — Ne, the general equation 

2 = AE— de (10) 
for nerve excitation is appropriate. Since E represents the stimulus inten- 
sity, E = o = am in our case, and with Aa = A, 


d 
qn Ain — 4 (11) 
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with A, and a as constants. Substituting (8) in (11) and integrating 
e= e “fey tM (Ri J+ ko) [a(1— co) — Re] }+Mke(kiJ + ko — @) (12) 
_ e—kIt+k,)t. Ma [RiJ =0g (RiJ + Ro) ] 
with 
AikiJ 


M = hy Eby eee 


where € is the value of € at = 0. 

For long exposures (¢ > ©) the afferent excitation attains an asymp- 
totic value 

€ = na "i be 

hi 


(13) 


where K = A,k2/a, which is a monotonically increasing asymptotic func- 
tion of J. Taking the initial conditions cp) = ¢) = 0 and denoting this by 
a circumflex (~) again, the afferent excitation for the att = 0 dark adapted 
system 


@=M{AikiJ [e-* — eGo th!) + he (kJ + ha— a) (1— e*) } (14) 


passes a maximum €» at 


1 kJ 


LCE) at Lely ey SORE ee a 


(15) 


whereas the greater J is, the greater the value of €,, and the smaller ¢,, 


will be. The upper part of Figure 4 shows the e, é-curves for two different 
light intensities. 
At ¢ = t, the condition for efferent excitation, 


e=h, (16) 


will be met. Then according to the one-factor theory of nerve excitation, 
the efferent excitatory factor e’ will be produced 


—=-= A’ (e—h) —a’e’. (17) 
Integrated, (17) yields 


¢ = eth Ai fo (e=2) ed rf (18) 
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with r = ¢ — ft and e'(t,) = ef. Upon quadrature and rearranging 
e = [Nei tr) + Pe hth) +r) 4.0, 


Mis (19) 
— eT [AN ee + P e~ hI +k) ty +O, — @/] 
1 
where 
yee 
NS Gg (OTM (hiJ + he) [a.(1 — co) — hel } , 
a Aja 
Sry heae oe Co (kiJ + Ro) | ; (20) 
Ay 
Q1= A thaM (hd + he 0) — A), 


and M as defined in (12). 


Neuron N; 


), 


Or 


Ficure 4. Time course of excitation in the afferent and efferent neurons, for the dark adapted 
system. Two different light intensities: Curves a and a’: After an afferent latent period 4 the 
é, ¢-curve intersects the threshold line h, starting efferent excitation (upper part). After an 
additional efferent reaction period 7» the 2’, /-curve intersects the threshold line h’, eliciting 
a reaction in the effector organ after a total reaction-time ¢ (lower part). Curve b: The @, 
t-curve does not reach the threshold h, thus no reaction can occur. Equations (14) and (19). 
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Equation (19) is the formulation of the time course of efferent excita- 
tion, the light intensity taken as constant, which acts on the effector organ 
as stimulus with a threshold h’. The similarity of the e’, écurve (Fig. 4, 
lower part) with the experimentally established action potential (Fig. 5) 
is obvious. It characterizes the sudden rising of the efferent excitation 
after an afferent latent period #, intersecting the threshold h’ of the effec- 
tor organ after an additional efferent reaction period 7» and thus eliciting 
a reaction after a total reaction time ¢ = ¢, + 79, the passing through a 
maximum of excitation, and the decreasing of the latter to an asymptotic 
value fori— o. 


PI 


Ficure 5. Time course of the action potential PIT after Granit and Riddell (1934) 


Results: The basic equations for a new theory of photoreception based 
on supposedly new, dynamical aspects have been formulated. The theo- 
retical time course of efferent excitation intensity agrees with experiment. 
So far no conclusions can be drawn as to the validity of this theory. In a 
forthcoming paper it will be shown that the formulation of reaction and 
latent times, intensity discrimination, and strength of reaction derived 
from the equations of this theory yields results fitting experimental data. 

The author is indebted to Professor W. Ludwig and the Zoological 
Institute of the University of Heidelberg for facilitating these studies. 
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It is shown that the validity of Branson’s integral description of metabolizing systems is 
subject to severe limitations. The validity is insured only in cases where the reaction is of first 
order, or quasi of first order. In all other cases Branson’s equation has to be modified to insure 
general applicability. The consequences of a different definition of the metabolizing function F 
have also been investigated. With the new definition F describes the pure effect of metaboliza- 
tion. It is found that in this case the integral equation is only capable of describing first-order 
reactions. With a slight modification of the integral equation it is possible to describe metabo- 
lites ‘‘with age,’’ which do not have reactions of definite order, but which satisfy the super- 
position principle. 

1. Introduction. In a series of papers H. Branson (1946, 1947a, b, 1948, 
1952a, b, c) has presented and applied a description of metabolizing sys- 
tems by means of an integral equation. A proper derivation of this equa- 
tion was never given, and neither were all the functions in the equation 
defined in an unambiguous way. The reader of Branson’s first paper 
(Branson, 1946) is furthermore confused by the first two examples given 
by Branson, in which the roles of accumulation and metabolization seem 
to be interchanged. It is claimed by Branson that the integral equation 
has a very general applicability. The main purpose of the present paper is 
to make a critical investigation of the validity of Branson’s equation, and 
it will be shown that there are severe restrictions on the applicability of 
this equation. 

For the following it will be necessary to refer to two of Branson’s papers, 
which we shall for brevity denote by B1 and B2 : B1 = Branson, 1946; 
B2 = Branson, 1952a. B1 is the first paper by Branson on the subject. 


The integral equation given in B1 is 
t 
M(t) =M(0)F() + f R(8)F (t— 8) a. (1) 
0 


Although this equation was never refuted by Branson, it was nevertheless 
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abandoned in B2, in favor of the following equation 
t 
M (1) = MoF (Mo, t) + _f R(8)F (M8, t— 0) d0 (2) 
0 


which differs from (1) in the argument of F. In equations (1) and (2) 
M(t) is the amount of metabolite at time t, M, is another notation for 
M(6), R is called by Branson the ‘accumulation function,’ and F the 
“metabolizing function.’”’ The functions R and F are not defined very 
precisely by Branson. However, it seems obvious, and has been affirmed 
by private correspondence between Branson and the author, that R(t) is 
the rate at which metabolite is entering the system from the outside, at 
time ¢. Whether or not R(é) is restricted to be non-negative is not discussed 
inigee de 

2. Two possible definitions of F. The definition of F is not so obvious. 
One is inclined to consider (2) in the special case R(¢) = 0 identically. 
The resulting equation: 


M (t) = MoF (My, t) 


describes the decrease of metabolite as a function of time as a result of 
metabolization, and the factor multiplying the original amount is exactly 
F(Mo, t). Thus, the metabolizing function seems to represent the pure 
effect of metabolization, i.e., the effect of metabolization if no additional 
amounts of metabolite enter the system. Underlying this reasoning is the 
assumption that for a certain metabolizing system F is a given function of 
its argument, or, in other words, that (2) should be valid for arbitrary 
M> and R(t) with the same function F. However, through a private cor- 
respondence between Branson and the author it became clear that F 
should not be considered as a fixed function. The functional form of F is 
allowed to depend on Mo and R(#). On the other hand, F has to satisfy (2) 
in a very special way: F(M,, t — 6) shall be the fraction, left at time ¢, of 
any amount of metabolite which entered the system at time 6. This gives 
at the same time the physical interpretation of F. It should be noted that 
in this interpretation the fate of individual molecules is followed. As a 
consequence, F could conceivably be determined experimentally by intro- 
ducing tagged molecules into the system (Branson, 1947a) and studying 
their disappearance in time. However, the existence of a function F with 
the required properties is by no means guaranteed, as will be made clear in 
Section 4. 

Recapitulating, we have two possible definitions of F , namely: Defini- 
tion (a): F is a fixed function of its argument and represents the pure 
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effect of metabolization; F shall satisfy (2) with arbitrary Mo and R(t). 
Definition (b): The functional form of F may vary with circumstances 
[given by My and R(é)]. However, F(M>, t — 6) shall have the property 
that it is the factor by which an amount entering at time @ is reduced at 
time ¢. For the sake of completeness we shall investigate the consequences 
of both definitions, but it should be kept in mind that only the treatment 
of definition (b) will be pertinent to Branson’s work. 

3. Consequences of definition (a). We will show that F(M,,¢ — @) cannot 
depend on M, if (2) is to be satisfied for all My and R(@) with the same 
function F. In order to show this we shall make use of the fact that in a 
physical system it should not make any difference whether a certain 
amount of metabolite is injected into the system just before or just after 
¢ = 0. Suppose an amount « of metabolite is injected into the system just 
before ¢ = 0, and no additional amounts enter after ¢ = 0. Then My = «x, 
R(@) = 0, and after substitution into (2) we have 


M(t) =«F (x,t). (3) 


Suppose now that an amount x is injected just after ¢ = 0 in the following 


way 
M, — 0 


NE Ce ey ee (4) 
T 
Rigi = 0 eos r 


As we let r — 0, this mode of injection approaches an instantaneous injec- 
tion of the amount x just after £ = 0. Substituting (4) into (2) we have: 


MW) =* [ FM, 1-6) 48, Pears) 
Fg 
M(t) =~ [ F Me, t— 0) 48, LE 
T Yo 
If we put 
t=ré, 6=77,; M(t) =%¢, M6@=yY,,; (7) 


we can write (5) as , 
ye=e fF (in 7 (E—1)) an. 
0 


Letting r — 0 we get: 
E 
j = % d . (8) 
lim Ve » [Fw 0) 7 
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The validity of this limiting process follows from the continuity of F and 
can easily be verified. From (2), by putting ¢ = 0, it follows that F(Mo, 0) 
= 1 for arbitrary Mo. Substituting F(y,, 0) = 1 into (8) we obtain 


lim ye= 8, (9) 


which expresses the linear increase of M from 0 to x when # increases from 
0 to 7, in the limit r > 0. 

In (6) we make the substitutions (7) for 0 and M,, and let 7 — 0. After 
taking the limit under the integral sign we may put y, = x7 by (9), and 
tn = 0. Thus we obtain 


. ; d t>0 
lim MW) =2 f Flan t) n, = 04 


and after putting x7 = y 
lim M (2) = [Fo, t) dy, t>0O. (10) 
Tm0 0 


We asserted that in a metabolizing system (3) and (10) should give the 
same result. Equating the two expressions for M(?) gives: 


oF (x,t) = fF (y, 1) dy, 1>0, 
0 


and differentiating both sides with respect to « 


F(x, 1) +2 PSB) pian, 

so that dF(x, t)/dx = 0 for arbitrary x. This proves that F(Mp, é) has to 
be independent of Mo, and F(M,, ¢ — @) independent of M, 9, aS was 
claimed in the beginning of this section. 

The first result we have obtained is that with the definition (a) for the 
metabolizing function equation (2) reduces to (1) in view of the fact that 
F(Mo,t — 6) is a function of t — 6 only. In order to obtain a further re- 
sult, we shall first make a distinction between metabolites whose metaboli- 
zation rate at time ¢ does not depend on their past history up to time 4, 
and metabolites whose fate does depend on their past history. An example 
of the case first mentioned is molecules undergoing a chemical reaction, an 
example of the second case is the decay of red blood cells, which have a 
more or less definite life span, like living organisms. For brevity, we shall 
speak of metabolites “without age” and “with age”’ respectively. If the 
metabolite is without age, the metabolizing function has to satisfy the cen- 
dition 

F(t) =F (6)F (t¢— 0) (11) 
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[see, e.g., B2, eq. (6)], for arbitrary @ < #. It follows from (11) that F(t) = 
exp [-ct], c a constant, which is the metabolizing function for a first-order 
reaction. Hence, if the metabolite is without age, and if we take definition 
(a) for the metabolizing function, then the integral equation (2) reduces to 
(1) and describes only first-order reactions. 

If the metabolite has age, then F(/) may have any form, as long as it is 
monotonic decreasing. However, it is no longer sufficient to know the 
amount Mo at ¢ = 0, since this amount could constitute a heterogeneous 
group of various ages. Equation (1) has to be modified slightly to the 
equation 


M(t) = f _R()F—0) 40. (12) 


In addition, R(¢) has to be non-negative. Whereas in the case of first-order 
reactions the fate of the metabolite can be described as easily by a differen- 
tial equation as by the integral equation (1), in the case of a metabolite 
with age equation (12) is the only possible description. 

4. Consequences of definition (b). The discussion in this section will be 
restricted to metabolites without age, which are the only metabolites con- 
sidered by Branson. 

In order to show the essential difficulty with equations (1) and (2) we 
introduce the function E(6, ¢). This function is defined as the fraction, left 
at time ¢, of an amount of substance having entered at time 6. With any 
given Mo, R(t) and metabolizing process, there will be associated such a 
function E(6, t). Equation (1) claims that E(6, t) = F(¢ — @), hence that 
E has a very special form in that it is only a function of t — 0. It will be 
shown that this is by no means true in general. Equation (2) claims that 
E(6, t) = F(M,, t — 4), or that E is a function of M, and ¢ — 0. This is 
not true in general either. If E(6, #) is not a function of ¢ — @ alone, it can 
always be considered as a function of @ and t — 6. Whether it is a function 
of M, and t — @ depends on whether 6 is a function of M,. The time 6 
will be a function of M;, if, and only if, M(t) is a monotonic function of /, 
either strictly increasing or strictly decreasing. Even if this condition is 
fulfilled, and therefore the integral description (2) possible, the usefulness 
of (2) seems doubtful. We shall return to this point later. 

An explicit expression for E(6, ¢) can be derived as follows. Let Q(¢) be 
the rate at which substance is being metabolized at time /, hence 


Q(t) = —~ EERO . 


Let x(9) be an (infinitesimal) amount of substance, introduced into the 
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system at time 0, and let «(#) be the amount which remains of it at time /. 
For convenience we shall speak of the substance which was introduced at 
time 6 as being tagged. At any time #, the ratio of tagged to untagged sub- 
stance is «(/)/M/(t). Under the assumption of homogeneous mixing and 
chemical indistinguishability of tagged and untagged substance, the rate 
at which the tagged substance is being metabolized is [x(t)/M (é)]O(2). 
Therefore, « has to satisfy the differential equation 

ax owe (t) 13 

a = ww Q(t) . (13) 
If we introduce the function 


(14) 


then the solution of (13) is 
« (t) =2% (6) exp [—S() +S(8)]. 
Since E(8, t) is by definition the factor multiplying «(6), we have 
E(6,¢) = exp [—S(@) +5(@)] (15) 


with S given by (14). 
That £ has to be of the form (15) with some S(é) follows also immedi- 
ately from the condition 


E (6)'t) =E (Ope) BW", 3); 0&2 SA 


which £ has to satisfy [analogous to eq. (6) in B2]. However, the explicit 
form (14) of S(¢) follows only after considering the differential equation 
for x(t). 

From (15) we can draw the following conclusion: In order that E be a 
function of ¢ — @ alone, it is necessary and sufficient that S be of the form 
S(t) = ct, with c constant. Consulting (14), this means that Q(é)/M(t) = 
c = constant. If M(é) is not constant, then QO(#) = cM(#) means that the 
substance is being metabolized according to a first-order reaction. If 
M(t) is constant, then it means that Q(¢) is also constant, and, as a result, 
that R(t) = Q(t) = constant. This is the case of a stationary state, in 
which the substance is supplied at the same rate as it is metabolized. The 
order of the reaction is then immaterial. In both cases the metabolizing 
function has the form F(t) = exp[—c#]. This form of F is natural in the 
case of a first-order reaction, but also in the case of a stationary state it is 
well known that a small amount of tracer disappears exponentially in 
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time. The differential equation for the tracer is of the same form as for a 
first-order reaction, so that we may consider the metabolizing reaction for 
the tracer to be quasi of first order. Hence, we have the result that only 
essentially in the case of first-order metabolizing reactions the function 
E(@, t) is a function of ¢ — 6 only, and has the specific form E(@, 1) = 
exp [—c(t — @)|= F(t — 0), so that F(t) = exp|—ct]. In all other cases 
E(6, t) cannot be a function of t — @ only. 

If O(¢)/M(?) is not constant, then E(@, f) cannot depend on ¢ — @ alone. 
It has been remarked already that E can then be considered as a function 
of 6 and ¢ — 6, and that F is a function of M, and t — @ if, and only if, 
M(t) is strictly monotonic. However, there does not seem to be a natural 
reason for wanting to replace a dependence on 6 by a dependence on M,, 
and it seems more natural to express the equation in terms of E instead, as 
follows 


M (t) = ME (0, t) + [RODE (S, i) dé. (16) 
0 


The function (6, #) certainly exists, even if R(t) is negative in some in- 
terval. Even though (16) looks like a more adequate description of the 
system than (2), its usefulness seems doubtful. One tracer experiment can 
only determine £ for one value of @. Furthermore, the form of the function 
E is not characteristic of the order of the metabolizing reaction. 

5. Recapitulation. In the foregoing it has been shown that Branson’s 
integral equation, with Branson’s definition of the metabolizing function 
[definition (b)], has no general validity. The validity is insured only in 
case the metabolizing reaction is truly of first order, or quasi of first order 
(for the tracer) in the case of a stationary state. In these cases equation 
(2) reduces to (1), the metabolizing function is a simple exponential in the 
time, and the system can also be described by a simple linear first-order 
differential equation. In all other cases the description by (2) is in general 
not possible. Instead, the general case may be described by (16), but its 
usefulness seems doubtful. 

If a different definition is made of the metabolizing function [definition 
(a)], such that F is a fixed function, then equation (2) reduces to (1) im- 
mediately, and describes only first-order reactions. If (1) is slightly modi- 
fied to equation (12), then F can be non-exponential, and (12) describes 
the metabolization of a metabolite “with age.” With the definition (a) of 
the metabolizing function F, this function describes the pure effect of 


metabolization. 
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The assumptions latent in the derivation of the integral equation of Branson are rendered 
explicit and discussed. It is shown that the equation is valid only for systems in which the 
substance disappears according to a linear rate law. 


1. Introduction. H. Branson (1946, 1947a, b, 1948, 1952a, b, c) has 
suggested and applied a mathematical description of metabolic systems 
based upon an integral equation. By a derivation analogous in detail to 
that of the equipment mortality equation (Churchill, 1944; Gardner and 
Barnes, 1942) the equation 


M(t) =M (0) F(t) + f-R(8)F (t— 8) d0 (1) 
0 


is obtained (Branson, 1946). Here M(0) and M(#) are the amounts of 
metabolite at times zero and #, R(t) is the “‘rate at which the metabolite is 
accumulating in the system,” and F(¢) is defined as the “metabolizing 
function.”’ In a later paper (Branson, 1952a) the functions F(t) and 
F(t — @) were replaced by F[M(0), ¢] and F[M(@), ¢ — 6] and although 
reasons for this change were not discussed, specific directions for con- 
structing these functions were given. It was also required that R(#) > 0. 

It is not clear that im general the reasoning underlying (1), or the 
modification just mentioned, is applicable to physicochemical systems, and 
the logical foundation of the Branson formulation has never been ade- 
quately discussed. In any event the advantages of the formulation, if they 
exist at all, have not been made obvious. Although there exists a brief 
discussion (Branson, 1948) of the relation between the Branson formula- 
tion and the usual method of differential equations, this relation has not 
been satisfactorily established. In fact the Branson formulation has been 
explicitly shown (Branson, 1946) to give known results for only one sys- 
tem, viz., the irreversible first order reaction. 

2. Purpose and scope of paper. Throughout the papers of Branson, 
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cited above, there is considerable lack of consistency with regard to the 
definitions and meaning of the functions F and R. Moreover, there is su- 
perimposed (Branson, 1952a, b) upon the basic formulation physical inter- 
pretation, which if ever correct is certainly not generally so, and the pre- 
scribed procedures for applying the formulation (Branson, 1947a, 1952a) 
are incompatible with what appear to be the current definitions of F 
and R. 

While it is not entirely possible and may not be desirable to separate the 
mathematical and physical aspects of the formulation, it is the purpose of 
this paper to discuss the logical and mathematical foundation of the basic 
formulation largely without reference to the above-mentioned physical 
interpretation and procedures which will be discussed elsewhere (Hearon, 
1953). In what follows it will be shown that the formulation is valid only 
when the chemical kinetic rate equation for the disappearance of the 
metabolite is 17(¢) = —gM(t) where g is required to be a constant. 

3. The mortality equation and the original Branson formulation (Branson, 
1946). Since the Branson formulation embodies essentially the same as- 
sumptions as does the equipment mortality equation, it is worth while to 
make clear these assumptions and their restricted validity for physico- 
chemical systems. 

In the equipment mortality equation, F(#) is the probability that a 
piece of equipment new (age zero) at time zero will survive to age ¢. The 
probability that a piece new at time @ will survive until TIME #, when its 
AGE will be ¢ — 0, is then F(¢ — 6). If M(#) is the number of pieces at 
time ¢ and new equipment is being put into service at the rate R(#), then 
M(t) is given by (1) in which the integral represents the number surviving 
at time ¢ from those put into service mew at times ranging from zero to the 
current time ¢. The function F(f) has the property that 


M(t) = (0) F(t) =n (0) F (t— 6) (2) 

where (7) is the number of mew pieces at time r. From (2) it follows that 
of F(t) 

M()).=M Oar ay (3) 


but it is not true in general that 
M(t) =M (0)F (i— @) (4) 
for, as can be shown directly from (4) or by comparison to (3), this implies 


that i 
F(t) =F (0)F (t— @) 2 . (3) 


which requires F(t) = e*', a constant. 
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Plainly, in the laws of transport and chemical kinetics there is no 
counterpart of age and the time of arrival of the molecules in the system.is 
a matter of indifference. The rates of flow and reaction are governed nut 
by the length of time which the molecules have resided in the system but 
by the instantaneous values of the gradients and concentrations. The 
original Branson formulation requires the validity of (4) and hence of Cs 
It is correct, therefore, only if the system is such that F(é) = e*. Accord- 
ing to the extended formulation, to be discussed in the next section, F(t) 
has this form only if the substance undergoes a first order reaction in the 
system. 

4. The extended Branson formulation (Branson, 1952a). In the extended 
Branson formulation, the function F(¢ — 6) is replaced by F[M(6), ¢ — 6]. 
It is shown that if F[4/(0), é] is constructed from the solution of the formal 
chemical kinetic rate equation for a reaction of definite order, then 


F(M (0),t] =F ([M(0), 6) F[M(6),t— 6], (6) 


which is the counterpart of (5), is actually satisfied. This modification cir- 
cumvents the difficulty entailed in (4) and (5). But this is not enough. 
There is also inherent in the Branson formulation this assumption: The 
system is such that of the amount, M(0), present at ¢ = 0 there remains 
at any time M(0) F[M(0), ¢] whether or not additional material has entered 
the system, so that if at the instant ¢ = 4 an additional amount M, is in- 
troduced there remains in the system at any ¢ > 4, the total amount 
M(0) F[M(0), ¢] + MM, F[M(h), t — th] and so on for any number of ad- 
ditions. The implication is that the initial amount undergoes reaction as 
if it alone were present in the system and each added amount undergoes 
reaction in a manner determined by the amount present at the time at 
which it was added but uninfluenced by subsequent additions. The total 
amount at any time, it is supposed, can thus always be resolved into the 
sum of amounts remaining from the initial and each added amount. This 
assumption can be shown to lead to definite and severe restrictions upon 
the metabolizing function and the rate equation from which it is obtained. 

Let the substance undergo reaction in the system according to the rate 
function f(M), i.e., in a closed system the substance, were it not formed 
from other substances, would disappear according to * 


AG OAG | (7) 


* This is in accordance with Branson’s procedure (1952a, p. 49). It appears later however 
(1952a, p. 58; 1952b, p. 60) that f is to include all negative terms or all factors, reaction and 
transport, contributing to disappearance. Only in a closed system, then, would f be the usual 
chemical kinetic rate law. This changes nothing in the argument to follow. 
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where c(t) is the amount (or concentration), the dot denotes time differen- 
tiation and f(c) < 0 for c > 0. The solution of (7) can be written as 


c(t) =G[c(0),t— 4] (8) 
provided that 
c(t) dx salar 
Joos Fay ks eae oe 


has a unique explicit solution for c(¢) and it can be shown that this is true 
if f(c) is bounded. The metabolizing function, by an obvious generalization 
of the Branson definition (Branson, 1952a, p. 49), is defined by the rela- 
tions 


0 0),#] =G([M(0O),t 
M (0) F [M (0), #] =G[M (0), #] \ (10) 


M(0)F(M(6),¢— 6] =G[M(@),t— 9]. 
The functions F and G have, by virtue of their definitions, the properties 
G(a,0)=a, F(a,0)=1, 
: } (11) 
G=f(G). 


The extended Branson formulation is then 


M(t) =M(0)F[M (0), 4] + f R(O)F (M8), t— 0] 40 


(12) 
‘R (6) 
M (6) 


The differential equation for M(t) is 
M (t) = f (M) +R() (13) 


G([M(6),t— #6] dé. 


=G(m (0), + f 


which expresses that the material disappears at the rate f(/) and appears 
or “accumulates” at the rate R(t) due to formation from other substances, 
flows into the system, or both. From (12), however, 


kame ‘R(6) = e 
M(t) =G[M(0), #] + Joareay G[M(6),t—6]d@+R(t). (14) 
Equations (13) and (14) agree only if 
f(t) = 1) + f FAX Git), 1-0} d0. (18) 


If R(¢)/M(t) is denoted by S(#) and the arguments of S and G are left 
understood, it follows, from (12) and (15), that (15) is valid only if the 
function f is such that 


i (G+ f'scde)=sG) + f'ss@ao. (16) 
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If f(M) = —kM, which corresponds to first order reaction, (16) is satis- 
fied. But under these circumstances G[M(0), ¢] = M(0) e-*, F(—#) is an 
integrating factor for (13) and the Branson equation is the formal integral 


of this linear equation. If f(MZ) = —k, which corresponds to zero order 
reaction, (16) can be satisfied only if 
¢ 
f SGd0=0 
0 


for all ¢. If this condition is satisfied the situation is trivial. 
Because of the relation between the functional forms of f and G, it is 
not an easy problem to determine what functions, if any, in addition to the 
linear homogeneous one, f(Z) = — kM, satisfy (16) for prescribed R(é). It 
can be shown however that for arbitrary R(#), only the linear homogeneous 
function satisfies (16). Let R(f) be such that at ¢ = #4 an amount M, is 
added to the system [i.e., R(#) = M,6(t — t,) where 6(¢ — 4) is the Dirac 
delta function]. Then (16) becomes 
fle (at (0), 1) +48 G14, tal b= {EL o), 41} 
(17) 
Ft Gla — a}, 


where A = G[M(0), t)]. In particular consider the time ¢ = h, then the 
equation to be satisfied is 


f(A +I) = f(A) +43 f(A). (18) 


Since M, is arbitrary it can be chosen as (nm — 1)A. Then (18) becomes 
f(nA) =nf(A), (19) 


which shows f to be homogeneous and of degree one. 

It cannot be denied that for a given functional form of f(/), a particular 
R(é) and initial condition can be found such that (16) is satisfied [e.g., 
the trivial case R(#) = M,5(t — t:), M(0) = O] but certainly for arbitrary 
f and R the Branson formulation is not valid (cf. concluding remarks). 

It is worth noting that, with proper modification, (16) is satisfied if 
f(M) = —g(t)M, where g(t) is any integrable function of ¢. It can be 
shown at once that if in (7) and (13), f(M) is replaced by g(t) f(M) then 
in (8), (9) and elsewhere ¢ — @ is replaced by 


h(t) —h(0) = fx (e) ae, 
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and the above development goes through unchanged terminating in (16) 
and its concomitant (19). In fact it is not necessary to assume that the 
variables are separable in (7), but this particular choice corresponds to the 
examples given by Branson and indeed to all known chemical rate laws. 
The above modification goes beyond the Branson formulation, for if g(¢) 
is not constant then F depends upon h(¢) — h(6) rather than upon ¢ — 8. 
In the Branson description it is insisted that F depends upon ¢ — 6 and 
therefore g(t) must be a constant. 

5. Concluding remarks. The crux of the Branson formulation and its 
defects is this: A function, denoted by F[M(6), ¢ — 6], is constructed such 
that if R(t) = 0 (i.e., no flow into the system and no formation from other 
substances) multiplication of this function by the amount present at time 
6 gives the amount present at any ¢ > 6. It is then asswmed that multipli- 
cation of this function [wherein M(6) is now the amount actually present 
in the system at ¢ = 6] by an amount added at time @ gives the amount 
remaining from this addition at any ¢ > 6. This assumption is correct only 
if the principal of superposition applies and this requires the linearity of 
the rate law or that F[M(6), ¢ — 6] not be a function of M(@) but take the 
form exp [—h(é) + h(8)]. For example consider the rate law (7) and let 
an amount M, be added to the system at ¢ = 4. Then according to the 


Branson formulation 
M,G(A,t-—t) 


M(t) =G(M(0),#] + 7A (20) 


where A = G[M(0), th] is the amount left at # = 4 from the initial amount. 
Now of the amount, A, left from the initial amount at ¢ = & the fraction 
G(A, t — ))/A is supposed left at any ¢ > 4. This same fraction applied 
to M, is assumed to give the amount of that addition which remains at 
t > t,. Of the initial amount there is supposedly left at any # the fraction 
G[M(0), t]/M(0). Therefore the initial amount and the added amount are 
assumed to disappear as if they alone were present in the system: If the 
amount in the system at time # decays to a given extent (determined by 
the assumed rate law) during the interval of time from time #; to the cur- 
rent time f, a portion of material added to the system at time f; is assumed 
to decay to that same extent during the same time interval. It is expected 
from chemical kinetics that this assumption will hold if the rate law is 
linear but not otherwise. The correct amount in the system at ft > # is 


which is obtained by straightfoward application of the rate law: There is 
present at ¢ = 4 the amount A, left from the initial amount, plus the 
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added amount Mj. At any subsequent time the amount is, according to the 
rate law, that given by (21). Equations (20) and (21) agree only if Gis such 
a function that 

G (Gla, t], d) 


GG [a,t] +b, ) =G(a, 
a ) (a,i+A) +56 Gia ; 


(22) 


is satisfied. This equation is satisfied if 
G(a,t) =a exp [—h(t) +h(0)] or if f(M) = —-A(t)M(t). 


It should be clearly noted that in the definitions (10), (6) is the actual 
amount at time @ present in the actual system for which R(é) ¥ 0. If it 
is true that R(é) = 0 then (6) is obeyed. For example for a zero order 
reaction 
K- (¢— @) 


F[M (6),#— 0). =1— (6) 


This function obeys (6) only if @(@) = M(0) — K®. In general (6) is 
obeyed only for a first order reaction. The definitions (10) are those dis- 
cussed by Branson (1952a). In fact in that paper it is asserted that if a 
substance disappears over several paths by reactions of different order 
then the F-function is a linear combination of the individual metabolizing 
functions the zero-order example of which is given above. There is an al- 
ternative point of view forwarded by Branson in personal communica- 
tions: F may be considered to depend, in its functional form, upon the 
particular kinetics of reaction avd upon the functional form of R(#). This 
alternative is discussed by Wijsman (1953). The advantage of such an 
alternative is doubtful. For then F is merely that function which for given 
M(t) and R(t) satisfies the Branson equation. The Branson equation then 
becomes a formula by means of which a function, the physical significance 
of which is not prescribed, can be computed. 

It has been pointed out in those cases for which the formulation is valid 
the Branson equation is simply the formal integral of a linear first order 
equation. In many cases then, the distinction between the Branson formu- 
lation and the method of differential equations does not exist. It is of 
course possible that there is, where R(é) is not known and only the form of 
F(é) is known, a distinct advantage in dealing with the formal integral. 
Most of the applications (Branson, 1952a, b) have been to systems for 
which it can be assumed that the differential equation is linear. These ap- 
plications have, however, been complicated by additional assumptions and 
questionable methodology which are discussed elsewhere. 
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BACTERIA* 
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A problem in probability is stated which includes the problem of the distribution of bacterial 
mutants as a special case. This problem is solved exactly but since the resulting expressions are 
too complicated for practical use, various approximate expressions for the distribution are con- 
sidered, especially for the bacterial mutation case. 


1.. It is well known that bacteria may become resistant to certain anti- 
biotics, bacterial viruses or other destructive agents due to a suitable 
“mutation” or a mutation of one of their ancestors. 

The problem of the statistical distribution of the number of mutants in 
cultures of initially non-resistant bacteria was first treated by S. E. Luria 
and M. Delbriick (1943) and since by D. E. Lea and C. A. Coulson (1949) 
and P. Armitage (1952). 

We treat a somewhat more general problem. The way in which this 
arises from the problem of the distribution of bacterial mutants is briefly 
sketched here; for more details see the author’s earlier paper (Tricomi, 
1952a). We consider a set of cultures of non-mutant bacteria which are 
allowed to grow for exactly m generations before the bactericidal agent 
is added. Now from each bacterium the following m + 1 results are pos- 
sible. 

1) It mutates in the first generation, giving 2”~! mutants in the final 
culture. 

2) Asecond generation descendant mutates, giving 2”? mutants in the 
final culture. 


m) A descendant in the last (mth) generation mutates, giving one 
(2™-™) mutant in the final culture. 

m + 1) The bacterium gives rise to no mutants. 

In general if 4, = number of mutants in the mth generation arising 

* Research sponsored by the Office of Naval Research while the author was at the Cali- 
fornia Institute of Technology. 
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from a mutation that occurred in the rth generation 
Dia for r= 1.2) occ 
and 
hm+i1 = 0, the case of no mutations . 

Now if p is the mutation rate, assumed constant, the probability of a 
mutation in the rth generation descendants of a single bacterium will be 
2 

We are here assuming that the mutations are so rare that we can with- 
out error count all the progeny among those subject to mutation. Then h 
is a random variable that can take on the values 


DD De oe ea AL 
with the probabilities, for the first m of these, 
p, 2p,..., 2p. 


Now let N be the number of bacteria in the original culture, where JN is 
of course very large, and put 


A= Z*™-I1N > - 


then our random variable / can take on (in addition to the value zero) the 
values h, forr = 1,2, ..., m with the corresponding probabilities \/N/A,. 
Our generalization consists in allowing the random variable / to take 
all the integral values 0, 1, 2,..., &, where & is some fixed integer, with 
the corresponding probabilities \,/. We are thus led to the following 
probability problem, which we state as an urn problem. 
We have a large number N of urns each of which contains many balls, with 


one of the numbers 0,1, 2,..., k on each ball. 

The probability of drawing a ball with the number h (h = 1, 2,..., k) 
from any urn is \,/N, where di, do, ..., ; are certain non-negative con- 
stants. 


We wish to determine the probability that when a ball is drawn from each 
urn the sum of the numbers on the balls will total n. In making this calcula- 
tion the Poisson probability distribution is to be used. 

In particular, the problem of bacterial mutation is a special instance with 
An = /h, where d is a constant and h has the range 0, 1, 2, 22,... te eos 

In this paper we give in Sections 2 and 3 an exact solution of this prob- 
lem, which was suggested to the author by M. Delbriick. 

Then, since this exact solution is useful for small values of only, we 
deduce from it two different approximate solutions: The first (Sections 
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4—6) under the hypothesis that & is small in comparison with n, and the 
second (Sections 7-10) in the bacterial mutation case where this hypothe- 
sis is not fulfilled but ) is generally large in comparison with 1. 

This last solution leads to some difficult but interesting questions of 
pure analysis, which are only sketched here. They are considered in more 
detail in another paper of the author’s (Tricomi, 1952b). 

2. Using the Poisson distribution the probability of drawing from the 


urns 7, balls with the number h(k = 1, 2,..., k) is 
eer 
r,! 


Consequently the probability that when one ball is drawn from each urn 
we obtain 7; balls with the number 1, 72 balls with the number 2,.. . , 7; 
balls with the number & is 


k 


k Th Th 
T[2he*= Sond ok At+Aet+.. - + Xx) 
Rae 3 Tre 


Th A=1 


and the required probability p, is given by the formula 


pr=eA>> TT, (1) 


where the sum is to be extended over all the sets of & non-negative integers 
11, 72, .-., 7% Such that 


mt2reot3rgt...ptkr,=n. (2) 


In particular we have 
Pome 4. 


On the other hand, if we consider the analytic function 


f(z) =exp(Aiz+ dos?+...+des*—A) 
we have 


k oo r oo k Th 

ie — an 

fea) mes Bien he 
where the last sum is also extended over all the sets of & non-negative 
integers 71, 72,..-, 7% Satisfying the previous conditions (2); hence the 
probability », coincides with the coefficient of 2” in the expansion of f(z) 
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in a power series. In other words these probabilities po, pi, f2, . . . have the 
generating function 
Dd) na" = f (2) =exp (zt do2?+.. . + AgsF— A). (3) 


3. This result has several important consequences. First from function 
theory we get the following two exact expressions for p, 


pa= =f (0) (4) 
and 
0+ 
ne ap te patel eA (5S) 


Secondly we can readily calculate the successive moments M,, of our 
distribution, beginning with 


My = Saas =1 (6) 


which is of course equal to 1. Then, since 


fet) =exp (pe! + ete ete Ay = > Sex 


n=0 


that is, 
t 
eee [Fare ) = 
we have ne: Leibniz’s formula 


5 (el) = a Be te aby (Are'+ 2r,e% +... + Rdg e*) | 


Sibi on (e*) = (Aref 2rrert hy, okt 
a eet apo e€ 16’ + 2rhe%*+.. + kd, ef) , 
and putting ¢ = 0 we get the recurrence formula 


Ma= (Gj) Mata ; (7) 


where 
An = M+ 22+... RN, Se has ht Mena Fs (8). 
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In particular we find in this way 
My = Ai, Mz = A2, + As, My = A8+ 3AyAo+ Ap , 
My = At+ 3A2+ 6AdA?+ 445A, + Ay... (9) 


If we take the mean M, = Aj as origin and call yu; = 0, me, ws, .. . the 
relative moments, the first few become simplified as follows: 


M2 = 0? = Ao, us = As, Ba = Agt 3.2, (10) 


Among other things this shows that Pearson’s parameters 6, and £ here 
assume the values 


U3 AS [ie ok 
ese ae we eee 
Sok ie Be ¥ xis (11) 


Similarly we can deduce from (4) a recurrence formula for p,. Putting 
A(z) =diz+ roz?+... + dr2*, (12) 


using the formula 
PF ay=F (a) A C27} 


and making m — 1 differentiations, we obtain 


400(2) = (G4) J) AOC), 


so that 
_ 2 ky fo (0) A) (0) 
=3DF ‘n—hl kh 
But 
f™ (0) s phgeimaiedl, 4 eam A) (0) me ee 
Ps Gr bie hi 0,(h+k), 
hence 
1 n* 
tn ==>) bMiPr-a (13) 
Nii 


where n* denotes the smaller of the two integers m and &. 
In particular we obtain in this way 


Po= e~4, pi=me 4, fro= (+5 d?) Cee at ae (14) 
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4. Despite the relative simplicity of the recurrence formula (13) it is 
not easy to obtain from it a general conception of our distribution. But 
this can be done by means of the method of the steepest descent applied to 
the integral (5) put into the form 


1 (o+) 
a Si 


er()dz,o(2) = A(z) — (n+1) log z—A, 
provided that the quantities defined in (17) below are all small in compari- 
son with 1. 

For this we have first to consider the equation 


1 
o'(#) = A'(s) ST Eno, 
which, putting 
Ai(2) = A’(2) = 2+2rAqs?+... +R S* (15) 
can be written 
A,(s) =n+1. (16) 


The path of integration must go through a root of this equation. Now 
considering that the integrand vanishes exponentially as z goes to infinity 
in the direction of the imaginary axis, provided that at least one of the 
\, with an even index does not vanish, we can use as the path of integra- 
tion a line parallel to the imaginary axis going through the unique* posi- 
tive root 2 of the algebraic equation (16). We obtain thus 


1 zt too 
em aaz J, HPL (a) +5; 6" (4) (2). 1 ds 


Qri 


Now we use the substitution 


= tiVertey 
and obtain 
a arV/2¢" (2) J “exp(—P+icsi+ cl . .)d é, 
where 
Signa o'"’ (£) tina g!Y (Zo) 
C= 3! ‘[é” (25) 1872" OR THC) TE (17) 


* The algebraic equation (16) has one and only one positive root fe because its left side goes 
monotonically from 0 to © on the positive real half-axis, 
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But 


Jew (-Ptiak + c+...) d= f e* tics 


—o 


— 3 
+ ct — cs +. .)di= Vx Gtia-Paf ae aly 
hence we have approximately 
p. wept AON Leos 
“ W2xdb"(z)’ 


provided that the quantities (17) are all small in comparison with 1. 
Now, since 
eb(%0) = @A(%)—A(l) g,—2-1 
we can write 
eA (2)—-A(1) 


Pa BVI Ag(h)” (18) 


where A(z) has the expression (12), 2) is the unique positive root of (16) 


and 
Az (Zo) = 22” (Zo). 


Moreover, since 


@'(s) = A"(2) —2E1 zo’ (2) = Ai(2) —n—1, 


zp” (2) +¢'(2z) =Aj,(2), 


we can also write 
Az ( 2) = 2A1( 2) = Zot 22r0Z,+. -- th dee. (19) 


5. The approximate formula (18) is still too complicated for most prac- 
tical purposes, but it can be substantially simplified in several cases. For 
instance, when 7 is large in comparison with & and with 1, it can be in- 
ferred from (16) that zo is also large, and of the order of »¥*, Moreover, 
since the quotients Ax(z)/A(z) and A2(z)/Ai(z) both approach & for 
z— o, if we put 


A (25) =F) Ao(%) =nke , 


it is to be expected that for large values of n, the quantities ’; and k, do 
not differ very much from the constant &. Similarly, if we put 


n \~77/*, 
es Ge. ; 
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that is, 


considering (16) it can be expected that LZ does not differ very much from 
the constant 1. 

These heuristic considerations aim only at showing that, when 1 is large 
in comparison with &, it is reasonable to expect that the probability #, can 
be approximately represented by an expression of the type 


e(n/ky)— Ear 
V/2rkon 43 : 


where fj, ko, L, and A are suitable constants. Now, if we put 2/k = v we 
have 


asia (zy el eA Tt 
Vr kon nh 24 WV kik V/ 2a ev” 


and, on the other hand, 
V2 ave ~v! =T (v1). 


Hence we are led to attempt to represent p, by means of an approximate 
formula of the type 


(20) 


where ki, K, and L are three suitable constants. 

6. Of course, all this has a practical utility if, and only if, it is possible 
to determine these three constants so that the first three moments of our 
distribution assume their correct values: 


M,=1, M,=A1, M2=Ai+As. 


Now this can be done, and in a very simple and elegant manner, pro- 
vided that a sum of the type 


YS n"p,, (were 051525 6. x) 


n=( 


can be identified with the (not exactly equivalent) sum 


= L wieee) 
Dd (4) "Ko = sD er 


v=0 
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In fact, since 


Ly Ss y ~ y 
ene Deter Deh era ere 


y=; ° v=1 7 v=1 v! 


we have the three equations 


Ke*=1, kKLe*=A,,  k&KL(1+L) =Aj+A,. 
But because of the first equation, the second can be written 
Ry — Ay ) 


and the third then assumes the form 
Ay (Ry + Ai) = Ai + Ag, 1k op Ry Ay = Ao. 

Hence 

Ae _M 


L 


k = — =e Ll 
7 Ai’ Ae K é ° 


We thus obtain the very simple approximate representation 
a3 ¥2 st withy > Ai 
pager, (v=tinL=%) (21) 
which differs only slightly from the classic formula of Poisson, to which it 
reduces for k = 1. 


7. The previous approximate solution is of no use in the particularly 
important bacterial mutation case, 


mas, ee aid De Ag ora ot d, = O otherwise , 


where } is a positive constant, because k = 2” now is generally a very 
large number. 


In this case, which requires special consideration, we have 
Ay(z) =d( a+ 2+...+ 2°") =F (2), (22) 
where F,,(z) denotes the mth partial sum of the well-known, non- 
continuable power series 
F(z) =2z+22+2'!+..., (23) 


which represents a function F(z) existing for lz | < 1 only. 
Such partial sums can be dominated using the entire function 


© 1 “5 } 
C1. eae & (24) 
n=1 


n!? 
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which the author has studied (including the case in which 2” is replaced 
by the power a” of any positive number a@ > 1) in the paper cited (Tri- 
comi, 1952b). By repeated application of the directly verifiable functional 
equation 
G(22) =G(z) + e—1, (25) 
we obtain the formula 
F,,(e) =m+G(2"5) —G(g). (26) 


Now, if we change ¢ into 2-™¢, we have 
F,(e® "*) =m+G(t) —G(2-!2). 


But G(z) is of the order of |z| since G(0) = 0, G’(0) = 1; hence as long as 
|¢| is bounded, we can assert that 


Fn (e-™) =m+G (1) +0(2-"), (27) 


8. In this case the approximate solution of the algebraic equation (16) 
for Zo is not easy. But we can abandon this and attempt instead to express 
both ” and p, by means of zo, or more exactly, of a parameter ¢ such that 


Zo = er ™t . 


In other words, instead of trying to express the ordinate », of the prob- 
ability curve as a function of the abscissa m, we try to determine a para- 
metric representation of this curve. 

For the expression of the abscissa the problem is very easy. From 
(16), (22), and (27) we obtain immediately 


n+1=dA[m+G (t) +0(2-*)], 
that is, 


n~d[m+G (t)]. (28) 


The expression of #, as a function of ¢ has to be deduced from (18) after 
the approximate evaluation of 


% d : 7 
A(m) —A(1) = f Ai(s) and Ay (2) = 2)Ai(%). 
But we have first 


Ay (go) =e? "Fi, (e? ™) = 2" 1G" (1) +0 (2-") ] 
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and, changing ¢ into e?-"", we find 


A (29) —A (1) raetf A,(e? ™4) du—2-» f- [m+G (uw) +0 (2-™) ] dx 
0 


= A2—™ [mi+G* (t) +0(2-™)], 


where 
t 
* f) = 7 . 
(1) = f G(u) du; (29) 
hence 
pw CAs _ exp [2-mi (— nt ym) + 02-"G* (i) $0 (2-™)] 
25 V 207A, (2,) V22d2™ [G’ (t) +0(2- m) | d 
that is, 
exp l= (%— Am) 2—- "i+ A2—-"G* (7) Ned 
oe 2/2 \/2r NG’ (1) (30) 
Now, if use is made of the change of coordinates: 
n=r(m+ 2x) x= m 
, (31) 
br= x9 Y = Nin 
and if we put 
2-™h=a, iG(t) —G*(t) =G,(é), (32) 


the reduced probability curve y = y(«) can be represented by the simple 
approximate parametric equations 


fae —aG,(t) 
S=iz), “= VE (33) 
which, in addition to #, contain the parameter a = 2-”) only. 

9. On the real axis G(t) is a constantly increasing function of ¢. More- 
over it is a completely monotone function, since all its derivatives are posi- 
tive. Furthermore the function G,(#) given by the second formula (32) is 
also always positive on the real axis. Figure 1 gives an idea of the behavior 
of G(é) for —8 < t < 6, while Table I gives some values of G(¢) as well as 
that of G’(#) and G*(#) in the range (—2, 2). 

As long as a is not too small the probability curve (33) can be easily 
drawn with the help of Table I. For instance, Figure 2, which shows the 
probability curve in the case a = 1, was drawn in this manner. 

However, if a is small (and in bacterial mutation case values of a of the 


20 


Ficure 1. The function G(?) 


TABLE I 
t POSITIVE b ¢ NEGATIVE 
|e 
G(t) G’(t) G*(t) G(t) G’(t) G¥(t) 
ieee be ky 0.000 00 1.000 00 | 0.000 00 0.000 00 1.000 00 0.000 00 
AN lie aie 101 69 034 06 005 06 | —0.098 36 0.966 37 004 94 
6s ea 206 86 069 61 020 45 193 52 935 60 019 55 
Ose Ann 315 67 106 74 046 56 285 62 906 14 043 53 
D4 earner 428 27 145 51 083 73 374 79 877 17 076 58 
0 Sohn 544 83 186 00 132 34 461 14 849 68 118 40 
ORG sae: 665 53 228 30 192 82 544 80 823 48 168 71 
7 Saeed aire 790 55 272 49 265 58 625 87 798 16 227 27 
Oh Sis ee te 0.920 09 318 67 331 07 704 47 773 74 293 80 
QO es te teers 1.054 35 366 92 449 76 780 68 750 57 368 08 
LEO cheers 193 55 417 36 562 12 854 61 728 21 449 88 


Ll el Lo 337 90 470 09 688 65 926 35 706 75 538 94 
Leyden de, 487 64 575 21 829 88 | —0.995 99 686 15 635 08 
USS Rs Nien 643.03 582 85} 0.986 16 | —1.063 61 666 37 738 07 
UE ere oe Bee 804 30 643 12 | 1.158 67 129 29 647 37 847 74 
85 iets oa Ete d OF E474: 706 16 347 42 193 11 629 13 | 0.963 86 
Te a 2.145 63 772 10 553 23 255 14 611 54} 1.086 29 
RETR 6, Got oe 326 26 841 09 | 1.776 76 315 45 594 76 264 84 
Udo c wemiiichs 513 96 913 26 | 2.018 73 374 11 578 57 349 34 
Oe Fig aug 709 03-| 1.988 78| 279 81 431 19 563 01 489 62 
PILE Fee ee 2.911 83 | 2.067 77 | 2.560 79 | —1.487 73} 0.547 99 1.635 53 


cen ne ee 
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order of 10~? occur) Table I gives us a small part of the probability curve 
in the neighborhood of « = 0 only. The remaining parts of the curve have 
to be drawn with the help of approximate expressions of G(#) both on the 
right (i.e., for £ > 2) and on the left (ie., for t < —2). 

The asymptotic behavior of G(#) is carefully studied in the paper cited 
(Tricomi, 1952b). From the results of this paper it follows in particular 
that on the right we have 


G(t) = e/2[1+0(e-“/4) ] G@ =} said eae (34) 


G* (t) =2e/2[1+0 (e-*/4) ] Gi (t) = e/? [t-—2+0 (e-*/4) 


FIGURE 2 


As a consequence, on the right there is no need to retain the parametric 
representation of the probability curve and we can represent it approx1- 
mately by means of the very simple equation 


—2axz(log z—1) : (x => 1) . (35) 


a 
pine V Tx : 
On the left the asymptotic behavior is less simple because one of the 
striking features of the entire function G(z) is the profound difference be- 
tween its behavior in the two half-planes Rez > 0 and Rez < 0. In par- 
ticular, on the negative real axis we have 


|t|G’ (4) =Q() +0 (4), (36) 
where the function 7(é) vanishes for |t| > © while Q(¢) is an oscillating 
function which practically does not differ from the constant 

A =log, e= 1.44270, 


because it can be shown that 
|O(t) —A| <0.0000143 . 
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Hence for t « —1 we have approximately (not asymptotically) 


G°@) ae G(t) ~ —loge|t] = — A log|t|, Gi () At] (37) 


and, consequently, the following approximate representation of the prob- 
ability curve 


t 
x= — A log |t|, oN ee (t< —1) (38) 


FIGURE 3 


holds. Or, if preferred, we can write 


Ab Eset =i) 
y Negi exp ( TA aAe . (+<0). (39) 


In case more precision is required, we can evaluate the function G(t) 
beyond the range (—2, 2) by means of the empirical formulae 


G(t) = — 1.995174 e+ e44+0.25/40.012, (422) 
bi (40) 
G(t) = — 1.33327 — A log (+) + 2ar100.25+02, (s—2) | 


where M = 0.43429. These formulae give generally satisfying results, i.e., 
errors of the order of 1%. 


10. Figure 3 represents the probability curve in the case a = 0.01. It 
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was drawn by means of the previous formulae (35) and (39) and of the 
numerical table for the neighborhood of x = 0. 

The main feature of this curve is the sharp maximum at about « = —5 
despite the fact that, according to the first formula (10), which now be- 
comes 

o7 =) (2™—1) ~ 22m, , (41) 
we have actually a large standard deviation. Hence a large standard devia- 
tion is not always the sign of a flat frequency curve. 

Another interesting feature of our curve, not without connection with 
the previous one, is the greatly decreased importance of the mean (x = 0) 
in comparison of the mode and the very pronounced skewness of the curve. 

Provided that a is so small that the mode ay lies in the range in which 
(39) is applicable as in the case a = 0.01, x can be readily determined by 
the condition 


Hi gettt  aae 
which gives 
%y= A log (2Aa) =1—m-+A log(AA). (42) 


The value of m corresponding to this value of xp is 
My = X+XA log (NA). (43) 


In the case in which m is so large that it is permissible simply to 
put m= o, the change of the parameter 2%» into ¢ is no longer suitable, 
and one has to use the formulae (16) and (18) where the polynomials A, 
A,, and A» must be replaced by the corresponding infinite series. 

It must be remembered that the validity of these formulas requires that 
(18) be valid. This requires that the quantities cz, cs, . . . defined by (17) 
be small in comparison with 1. Now, in the bacterial mutation case, we 
have 

ée=0 (i), 6Ge= OK 2), 2 223 
hence \ must be sufficiently large, perhaps of the order of 10. 

When this condition is not fulfilled, we must return to the recurrence 

formula (13), which in the bacterial mutation case becomes 


n= DD, Pa—2', = v=min(2™"}, [log, ~] ), (44) 


where (log, ) denotes the greatest integer in the logs 2. In other words, 
we have to stop the sum at such an hy that 2% < m but 2hot! > n, as long 
as this fp is less than m and at m — 1 in the contrary case. 
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In this way we find in particular, supposing m = 4, 


1 
bo=exp[—2d(1—2-™)],  pi=Apo, p2=zZr(1+2) hos 
P= ZMOF3) Py, P= GAO +ON+INES) fo, 

1 


Ps MM 10d? + 15+ 30) 95 ee 


For the sake of brevity we omit the proof of some further properties of 
our distribution, for instance, the fact that in the case k = 2 the probabil- 
ity p, can be exactly expressed by means of an Hermite polynomial of 
imaginary argument. 
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Simple reaction and discrimination reaction, under the influence of imitation, are considered 
for the situation in which the stimulus or the stimuli vary slowly with time. The result is 
analogous to hysteresis under certain conditions. The calculations are facilitated by the 


solution of 
at+px 
z= f- g(t) dé, 


g(&) being the normal error function. Values of x(a, 8) are given in a table. 


Simple reaction with the effect of imitation. The simplest structure that 
one can consider to represent an individual is that of a pair of neural ele- 
ments in which a stimulus S; acts upon the first element, while the second 
element acts upon some effector element leading to a reaction Ry. In 
addition, since we are interested in the case in which the behavior of one 
individual may affect that of another, we shall suppose that the reaction 
of each individual is a stimulus affecting the behavior of the others. Fur- 
thermore, we shall suppose that the individuals have identical parameters, 
except perhaps for threshold values. For convenience, we shall consider 
that the thresholds are the same also but that various factors tend to add 
to or subtract from this value, and that these effects are randomly dis- 
tributed. 

Let ¢ be the excitatory factor between the pair of neural elements repre- 
senting an individual, e being produced by the first and acting upon the 
second neural element. Let # be the threshold of the latter. Let ¢ be a 
measure of the effect of random variations with standard deviation ¢, 
which subtracts from the effect of ¢. Since we are interested in situations 
in which the behavior of one individual is influenced by the reactions of 
others, we must suppose that somehow the value of ¢ is changed as a 
result of the behavior of others. Let the amount added to € be proportional 
‘to the number WN of other individuals which show the reaction. Then if } 
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is the coefficient of proportionality, we have (Rashevsky, 1948, 1949; 
Landau, 1950; Rapoport, 1952) 


ie ABS) ene (1) 


the quantities A, Z, and a being the same as those defined by N. Rashev- 
sky (1948). If S varies slowly enough, 1.e., dS/dt « aS, then to a satisfac- 
tory degree of approximation the following relation is true 


AA ¢ ON 
0, ps Bs) ae ae (2) 


If we consider all individuals to be equivalent, then the number W divided 
by the total number Wr is the fraction exhibiting the response. This frac- 
tion will be very nearly equal to the probability P; of the reaction. If 
g(€) is the Gaussian error function according to which ¢/ac is distributed, 


and if 


Gi= fede, gl) =e, (3) 


then P; is given by 


jer = 


[A E(S,;)+bN—ah] /oa 
ie g(£)dE=346 


(2e + oN *), (4) 


od 


—o 


h and o having been defined above. The range of £ in expression (4) is 
obtained by setting de/dt = 0 and requiring that ¢ be greater than h. If 
we set P; = N/Nyz, and if 


aly FESO) _& (5) 
ac Co 
and 
_bNy 
Rela meet (6) 


then we find from equation (4) 
P;=3+G [a (t) + 6P;). (7) 
From an inspection of (7) and (3) we see that we cannot solve for P; 
directly. We shall find it convenient to define a function H (x, y) as follows: 


a+yH (x, y) 
W(x, 9) =f g(t) de=b4G6 let yH (a, y)]1. (8) 


By comparing (7) and (8) we see that P; is given by 
P,(t) =H [a(t), B). (9) 
In Table I values of H(x, y) are given. Note that for |y| less than about 
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1 H(x, y) is satisfactorily approximated by 


fsa 
be ye) 


10 
H(“,y) = a 


[1—axste (x) fee) ae]. 

Now suppose that S increases slowly with time until a time ¢ = T, that 
is, for some reason the stimulus which tends to produce the response be- 
comes progressively larger for a period of time. We may follow the per- 
centage reaction as a function of time by means of (9), with a(t) being 
given by (4). For simplicity, let E[S(é)] be linear in the time ?#. If B < / I; 
then it can be seen from expression (7) and its derivative that it has only 
one root so that P; is single-valued. Thus we obtain for P;(¢) a relation of 
the kind illustrated by the solid curve in Figure 1. For this curve 6 


1.00 


.50 


Ficure 1 


= 3, h/o = 2.25, r = AE(S;)t/ac, being a measure of time and £ being 
the derivative of £ with respect to time. The dotted and broken curves 
will be discussed subsequently. 

If B > \/2zn, there may be three roots of (7). This situation will occur 
if 8 = —2a and B > »/2z. There will also be three roots when 8 is suf- 
ficiently close to — 20 (cf. Table I). It can be shown that the middle root is 
an unstable point. Two examples in which 8 > 1/27 are shown in Figure 
2. In curve I the proportion of individuals exhibiting the reaction increases 
slowly, then at r = .82 it jumps discontinuously to a high value, after 
which it rises slowly. Now if the stimulus should somehow be reduced 
slowly, the proportion would at first fall slowly. Then at r = .67 it would 
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fall abruptly to a low value, returning slowly to its starting value. This 
case was obtained by setting 8 = 3, h/o = 2.25. 

In curve II of Figure 2 imitation results in a permanent change in the 
behavior pattern after the removal of the external stimulus. In this case 
the parameters are the same as in curve I, except that B = 5. 

Discrimination reactions with one reaction affected by imitation. Consider 
the situation in which there is a choice between two reactions R and Ry. 
Let S and S; be the stimuli which normally lead to R and R; respectively. 
Let the reaction to R; be influenced by imitation such that the strength of 
the effect is directly proportional to the number of individuals exhibiting 


1.00 
0 
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4 
oe, 
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the response R;. If S and S; lead to a simple discrimination mechanism 
with threshold h’ (Landahl, 1938, 1950; Rashevsky, 1949; Landau, 1950) 
we find for the probability of responses R and R, the following expres- 
sions: 


Py=34G [a (0) +6P:—') (11) 
d 
a P=$46[—a'(t) —6P,—#'), (12) 
with : 
dey 
ef = A2 is, y1-“ts@1, WAG (13) 


and 6 again given by (6). : re 

If a’(é) is, for example, linear in time and changing slowly, we fin i 
given by expression (9) in which a(f) is replaced by a’(t) — h’. Then P is 
given directly by expression (12). 
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In Figure 1 the curves may also be used to represent this case if 
h’ = 0.25, 8 = 4,7 = AE(S)t/ac, AE(S)/ao = 2,80 that o’(t) = —2+7. 
The solid curve represents P; while the broken curve represents P. The 
dotted curve represents the probability Py of no response. Similarly the 
solid curve of Figure 2 represents P; for the same choice of values of h’ 
and a’(t), but with 8 = 3 in curve I and 6 = 5S in curve II. The corre- 
sponding values of P and Py are not shown in this case since the figure 
would become overly complicated. 

Another situation to which these results may apply is that in which S 
remains constant but the population density changes with time. If, over a 
period of time, NV increases approximately linearly, then we obtain results 
similar to those shown in the figures. If a is too large, the variation of the 
fraction exhibiting the response will increase continuously, but if a is 
small enough, there will be a discontinuous change and the fraction at any 
time will depend upon whether WN has been increasing or decreasing and 
upon the past values of JV. 

Discrimination reaction with both reactions influenced by imitation. If 
both reactions, now denoted by R; and Re, are influenced in a similar way 
by imitation, it can easily be seen that we have the following expressions 
for P; and P, 

P,=34+G [a’ (t) +i Pi— B2P2), (14) 


P,=3+G[-—a’ (t) —2h’ —B:Pi+ B2Po) , (15) 


8, and B, being given by expression (6) with b equal to 2 or 62 in order to-_ 
allow for the possibility that the coefficients of imitation may be unequal. 
If h’ = 0 so that there is no threshold involved, i.e., one of the two reac- 
tions always occurs, then P; is given by 


P, =H [a’ (t) — Bo, Bit+ Be), (h'=0) (16) 
with Pp = 1 — Py. 

If h’ # 0, P, and P» for any given ¢ may be obtained by the intersection 
of two curves. The first curve is the result of plotting the values of 
H{a'(t) — BePs, 6] as the ordinate against P. as the abscissa. The second 
is the result of plotting the values of H[—a’ — 8,P,, 8.] as the abscissa 
against P; as the ordinate. Since the curves are not necessarily continuous, 
extra care must be taken in the interpolation required to find an inter- 
section. 

For the case in which h’ = 0 we may use only the difference in the num- 
bers of individuals showing reactions R; and Rp. Such a situation has been 
discussed by N. Rashevsky (1949) and H. G. Landau (1950). If X = NrP; 
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and Y = N;P2 are these respective numbers, then from (16) we find 
X—Y=2NfH [a’ (t) — By, Bi + Bo] —Nr. C17) 


If a’ is a constant and 8, = fe, the above expression, together with ‘Table 
I, gives the difference between numbers in the two groups of individuals 
after equilibrium has been reached. 

In view of (16) it can be seen that 


Pee a Ce, (18) 


This equation enables us to obtain values not listed in the table. It may 
also be used to obtain values otherwise requiring interpolation. From it 
can be seen that when x = —x% — yor y = — 2x the value of the function 
and its complement are both roots of equation (8). When y < 1/2z7, this 
root and its complement are 3. For y > 1/27, there are, in addition to the 
root 3, two roots whose sum is unity. 

It has been assumed as a first approximation that the effect of imitation 
is proportional to the number JW of individuals who show the reaction. 
Actually one would expect a more complex relation. In terms of our model 
the influence of the behavior of other individuals upon that of a given one 
can be determined as follows. Let f(m, 2) be the amount added to ¢ due to 
the observation of the individual of m “‘correct” and  “‘wrong”’ responses. 
For simplicity we shall assume that each response observed is that of a 
separate individual. If Ac is the difference, e, — ¢€,, between the excitatory 
factors due to the “‘correct’’ and the ‘‘wrong”’ responses in the absence of 
other individuals, and if / is the threshold used above, then the probabili- 
ties of “correct,” “wrong,” and “equal” responses after the individual has 
observed m “correct” responses and n “wrong” responses are respectively 
(cf. Landahl, 1939) 


pieagte (Atlin —*), (19) 
peag—o (tin a') (20) 
Fret cyl bie! sie 4 C29) 


Here the prime denotes the three-category case; the asterisk denotes the 

: * 
probability as altered by imitation. In the two-category case, i.e., P7 = 9, 
the probabilities are determined from 


pent—ps=3te (tL). (22) 
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If P, and P, are the corresponding values of P and P*, in the absence 
of other individuals so that f(m, 2) = 0, if G(y) is defined by 


GUY) = Rtn ol = GCE) S (23) 
then, since Ae/o = G"(P. — $), 


£5 (m,n) =G-\(P¥— 4) —G(P,-4), (24) 

* 5 (m,n) =G7(P*— 4) GP 4), (28) 
and, also, 

* 5 (m,n) =G-\(§ —Pi*) =G71(4-P)) (26) 


If the latter values calculated from experimental data do not agree with 
(24), then /# has been altered by imitation. If various values of Ae are 
used, the effect of Ae on f(m, m) can be determined. Note that the three 
values of G/o are experimentally independent. 
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discussing the manuscript. The author is indebted to Dr. H. G. Landau 
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On the basis of simple physical considerations the blood flow in a branching circulatory 
system is studied. The case of two groups of parallel vessels is treated. The vessels of the same 
group are supposed to be identical. The resistance of each group is determined by the resistance 
of each vessel in the group and by the number of vessels in the group. From the dependence of 
the resistance of each vessel on its radius an expression is obtained for the blood flow through 
each group of vessels in terms of the numbers and sizes of the vessels in each group. The num- 
ber of open vessels in an organ and the radius of each of those vessels are assumed to depend on 
the metabolic rate of that organ. The relations so obtained, together with the expression above, 
are applied to derive the blood flow through an organ as a function of the metabolic rate of that 
organ. It is indicated that the relations obtained might describe the shifting of blood from one 
organ to another if the activity of one of them changes. A way is pointed out to treat neural 
regulation of this phenomenon. 


The branching of vessels, in particular the smaller ones, e.g., the capil- 
laries, gives rise to rather complicated problems. Because it is a character- 
istic feature of the vascular tree it is of great importance in any theory of 
blood circulation. This structural aspect of the circulatory system is, for 
example, very important in the determination of the blood flow in the 
several parts of the system. This is not only true for the resting state, but 
also if changes in the activity of one or some of the organs occur, or in 
general under conditions of local need of higher blood flow. This shifting 
of blood in the lower extremities has been observed by H. L. Reese eé¢ al. 
(1952), who found that an increase of the skin blood flow produced, e.g., 
by heat, or a vasodilator such as priscoline, is followed by reduced calf 
muscle circulation. They observed also that when the circulation in the 
anterior leg muscles is increased by exercise, the circulation in the an- 
tagonistic calf muscles is reduced. They also pointed out the clinical impli- 
cations of these studies. M. Friedlander et al. (1940) found that muscle 
temperatures usually did not parallel changes in the skin temperatures. 
They concluded that in a given extremity the skin circulation and the 
muscle circulation are under separate control. O. Horwitz et al. (1949) 
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observed no significant change in the cardiac output after the use of heat 
or priscoline. The shifting of cutaneous blood flow has recently been em- 
phasized by M. E. De Bakey et al. (1947). They have called this phenome- 
non the “borrowing and lending phenomenon” or hemometakinesia. 

As a starting point for the mathematical treatment of such problems 
we will first study the simple branching system sketched in Figure 1. If 
we denote the pressure at the branch points by P’ and P”, the resistances 


«i ee 
POs yeh, 
I— P’ p”  I— 
edz iat igemehegiy 
FIGurRE 1 


of the branches by R; and Re, the mass flows in the branches by J; and 
Iz, and the total mass flow by J, we have 


P’—P"=R,i,, (1) 
P'’—P” =R.I2, (2) 
whereas the continuity of mass flow gives 
T=h+le. (3) 
From (1) and (2) we find 
I, R 
jay (4) 
Solution of equations (3) and (4) for J; and J» yields 
he ean 7 (5) 
i= RE i (6) 


For the case of Poiseuille flow we have 


R, = — and Ro = ——— 
T 


in which denotes the viscosity of the fluid, / the length of the vessels, and 
r, and ro the radii of the vessels as indicated in Figure 1. For the case of 


BLOOD FLOW 303 


vessels for which Poiseuille’s law does not hold we will assume 


G Gar 
ae and R= with a>1. (7) 


1 2 


The special case in which a = 4 and 


eqslih 


Tv 


CG 


corresponds to that of Poiseuille flow. 
Introducing (7) into (5) and (6) we obtain, after simple rearrangements, 


re 
Ents . (9) 


Introducing the equivalent resistance R, of the parallel circuit defined by 


P! —P" =R,I (10) 
we easily derive from (1), (2), (3), and (10) 

epee 

yA PD ee ay 


If both branches have the same resistance R we find from (11) 
ay 


aoe (12) 


and, in general, if there are ” parallel branches of the same resistance R 


Let 
eR (13) 


Let us apply this result to the case in which we have two groups of parallel 
vessels. The first group consists of m vessels of radius 7, the second group 
of n vessels of radius 72 (Fig. 2). The two groups may represent the vessels 


in two organs. 
- According to (13) the total resistance Rj, of the first group of m vessels 


each with resistance R, is given by 


Ry =— (14) 
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and that of the second group of vessels each with resistance R, similarly 


by 
R 
Roy =~. (15) 


The total blood flow J, in the first organ is given by (5) if we substitute 
R,, for R,; and Re, for Re 


Ro 


ieee (16) 
Sie pm psig Fs 
Pp’ p” 
Ficure 2 
Similarly for the second organ 
__ Rw 
ln pea gets (17) 


Elimination of R; and R» from equations (7), (14), and (15) and substitu- 
tion of the expressions so obtained for Ri, and Ro, into equations (16) and 
(17) yields, after some simple rearrangements, 


_ mrt 

are are ; (18) 
_ aes 

tt ara i (19) 
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It is easy to extend these formulae so that they correspond to more 
complicated cases than the one illustrated in Figure 2, e.g., the case in 
which the groups of parallel vessels are fed by two vessels which are them- 
selves branches of the main vessel in Figure 2. In such a case the pressure 
differences of the two groups will, in general, be different. However, for 
the sake of simplicity we will restrict ourselves to the case of Figure 2. 

In general 7; and will be functions of the metabolic rate gi of the first 
organ, whereas 72 and m will be similarly dependent on the metabolic rate 
gz of the second organ (cf. D. L. Cohn, 1953). We will assume here that 
these relations are given approximately by 


11 = Tot 4 (Gi — 0) 5 (20) 
12 = Tay + a2 (g2— goo) (21) 
m = m+ b1(g1— Gio) 5 (22) 
N = Na+ be (g2— Gao) - (23) 
The quantities which have an index 0 refer to the resting state of the 


organs. 

In general we suppose that there exist g’s which are related to metabolic 
activity, such that the relations (20) through (23) hold; they may be color 
changes of metabolic indicators. 

Even though the pressure difference P’ — P” will in general change if 
the metabolic activity of one of the two organs above changes, P’ — P”” 
will remain the same for both organs in a structure as represented in Fig- 
ure 2. The latter will be the case if the blood flow to both organs is supplied 
by the same vessel. In this connection we would like to point out that 
according to the results mentioned above of Horwitz et al. (loc. cit.) the 
pressure difference P’ — P’’ probably does not change if the skin circula- 
tion is increased by heat or priscoline. 

To simplify further notation we will introduce the following abbrevia- 


tions 
gi — Gu =A (24) 
and 


Jz — Yoo = AQ - (25) 


Substituting (24) and (25) into (20)—(23) and the results into (18) and 
(19) we obtain 


os, (mip t+ b:Aq1) (rio + LN eee | we 0) 
1 “Cryo + 61091) (710+ 41491) *+ (M29 + 2A go) Cant aq)" 
(129+ beAge) (120+ a2A gz) * EpGhi) 


as (mip t+ 01091) (riot G1A91)*+ (M20+ boA q2) (rap + @2A qa) * 
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In case the first organ is the skin, we may suppose that the resistance 
changes of that organ are primarily due to changes in radius of its vessels, 
whereas the number of the vessels does not change much (cf. M. van 
Dobben-Broekema and M.N. J. Dirken, 1950). Therefore we will assume 
b; = 0. On the other hand, if the second organ is the muscle, we may sup- 
pose that the resistance of that organ primarily changes due to alterations 
in the number of the vessels (cf. Krogh, 1929, p. 63), whereas we will 
consider the radius of the vessels to be constant. This leads us to also 
assume d2 = 0. 
The blood flow J, through the skin is then obtained from (26) with 
Oo = & = 0. 
T= Mo (Tio 414 41) * I 
* mig (rio + aA g1) *+ (Moy + boA qo) fre, (28) 


and similarly (27) with b, = a, = 0 gives the blood flow J,, through the 
muscle 


i (M9 + boA go) 729 I 
"mao (rip t Agi) *+ (29+ boA Qo) ae (29) 


For the case of exercise of the muscle without increase in skin metabolism 

(28) and (29) give, with Ag; = 0 and Ag, = Agn, 

a. = are 
i-y(s- pAg ia 


a BAG) 
"“T+7A+BAgn) -’ 


F (30) 
(31) 
with 

p=— (32) 
and 

yaa zn. (33) 


Mio \ Tio 
The plot of 7, and J,, given by (30) and (31) as a function of Agm, the 
increase in the rate of the metabolism of the muscle, is sketched in Figure 
3. This figure is based on the linear approximations assumed in (20)-(23). 
However, in general, J will also increase if gm increases. If this is the case 
we easily derive from (30) and (31) 


Im 
J ~ VU +BAga) (34) 


According to (34) the plot of 7,,/Z, against Agn is a straight line (Fig. 4) 


ow, Re ea) Whisk, GO Lee 
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within the range of validity of the linear approximation (23) and the 
assumptions made. 

For the case of a change of the skin metabolism, as, for example, might 
occur with change in the skin temperature, without change in the muscle 
metabolism, we have Ag; = Ag, and Ag. = 0 and substitution of these 
values into (28) and (29) yields 


(1+ dAq,)* 


I,= 5 
(+ oaq)*+y-’ S22 
eae ae wh es 

Se 6G TCR ae sae 
with y given by (33) and 
a (37) 
Tio 


This dependence of J, and J,, given by (35) and (36) as functions of Ag, is 
sketched in Figure 4, a being >1. This figure is also based on the linear 
approximations assumed in (20)—(23). 


O AG on 


Ficure 3. Blood flow Im through muscle and J, through skin as function of the increase 
Agm of metabolic rate of the muscle as given by (30) and (31). 


AS mn 


FIGURE 4 
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The curves in Figure 3 and Figure 4 are determined by the two parame- 
ters B and y. The latter can be determined from the intercept and slope of 
the plot of J,,/I, against Agqn. After having determined vy the curves in 
Figure 5 depend on the two parameters a and 6. These parameters can be 
determined, at least in principle, from the plot of either J,, or J, as a func- 
tion of Ag,. This would provide another method to obtain the exponent a 
of the power of the radius in the expression for the resistance, as is seen 
from (7). The determination of the parameter y would give the ratio of 
the effective resistances of the two organs in the resting state, as is seen 
from (33). The value of the parameter 6 would give a measure of the rate 
of opening of capillaries in the muscle with insrease in metabolic rate of 
that tissue. Similarly the value of 6 would give a measure of the ability of 


O 
Ags 


Ficure 5. Blood flow Z through muscle and J, through skin as function of the increase 
Aqs of metabolic rate of the skin as given by (35) and (36). 


the capillaries of the skin to dilate upon increase of metabolic rate of that 
tissue. Taking for the g the rate go, of oxygen consumption, the above rela- 
tions might give the dependence of the distribution of the blood flow on 
the oxygen consumption. Neural stimuli might change a; and by. 
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The partial differential equation of the random walk problem with persistence of direction 
and external bias is derived. By persistence of direction or internal bias we mean that the prob- 
ability a particle will travel in a given direction need not be the same for all directions, but 
depends solely upon the particle’s previous direction of motion. The external bias arises from 
an anisotropy of the medium or an external force on the particle. The problem is treated by 
considering that the net displacement of a particle arises from two factors, namely, that neither 
the probability of the particle traveling in any direction after turning nor the distance the 
particle travels in a given direction need be the same for all directions. A modified Fokker- 
Planck equation is first obtained using the assumptions that the particles have a distribution of 
travel times and speeds and that the average time of travel between turns need not be zero. 
The final equation incorporating the assumption of a persistence of direction and an external 
bias is then derived. Applications to the study of diffusion and to long-chain polymers are then 
made. 


Introduction. The classical random walk problem deals with a particle 
moving in a series of steps, where the length of the step, time between 
steps, and direction of the step are independent of each other and of pre- 
ceding steps. This problem was first formulated by K. Pearson (1900), 
although it had been treated and solved in a different form as early as 1880 
by Lord Rayleigh (1945). For some excellent reviews of the history and 
work done on this problem the articles by S. Chandrasekhar (1943), G. E. 
Uhlenbeck and L. S. Ornstein (1930), L. Infeld (1940), M. C. Wang and 
G. E. Uhlenbeck (1945), and W. Feller (1949, 1950) may be consulted. 
Classical random walk, as defined above, may be considered to be a Mark- 
off process, and as such may be described by the Fokker-Planck equation 
(Kolmogoroff, 1931), the solutions of which yield the probability of a 
particle being at any point at time ?. 

Random walk has been applied to numerous different physical situa- 
tions, many of which are of interest to biology. In the diffusion of gases, 
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a method analogous to the classical random walk treatment was used and 
was thought to yield the correct results (Meyer, 1899), but it has been 
shown that this method is not really applicable (Chapman and Cowling, 
1939; Jeans, 1940; Furry, 1948), since the postulated independence be- 
tween the steps is not valid for gas particles. However for diffusion in 
dilute solutions and for Brownian motion, the random walk method does 
seem to be applicable and yields useful results (Einstein, 1926; Jost, 1952). 
In the treatment of genetic variability in populations (Wright, 1945; 
Feller, 1951), the random walk method has been found to be applicable. 
However, as yet the equations have not been capable of solution for any 
but a few restricted cases. 

Another interesting application of the random walk method is found in 
the study of long-chain molecules. One may perhaps assume as a first ap- 
proximation that there is no fixed bond angle between atoms. This as- 
sumption is of course false, but one may consider a group of atoms as a 
unit, and this assumption will hold fairly well for the unit. Then by con- 
sidering the chain as having been formed by a particle executing a random 
walk, where the atoms of the molecule are to be placed at each point where 
the particle made a turn, the probability distribution of the possible 
lengths of the molecule may be calculated (Guth and Mark, 1934). Thus 
upon stretching the molecule (as in rubber or muscle), the entropy changes 
may be calculated (Treloar, 1942). 

As a final example, it may be pointed out that the motion of organisms 
under the influence of various stimuli (Fraenkel and Gunn, 1940) may, to 
a good approximation, be treated by this method. This has been done by 
K. Pearson (1906) for the case of random migration of animals without 
any stimulus, and the results were applied to the migration of mosquitoes 
into a new region. A more general treatment has been given by J. G. 
Skellam (1951), who took the population growth into account. The ran- 
dom walk method has also been used by C. A. Coulson (1947) in an appli- 
cation to the movement of some Mollusca in the presence of the stimulus 
of non-isotropic lighting, and by D. H. Wilkinson (1952) in a discussion of 
a theory of pigeon-homing by means of random search. 

A more generalized problem may be considered, however. That is, the 
individual particle’s steps need not be independent of each other, i.e., the 
random walk is not a Markoff process. The most familiar example of this, 
although applied to a different situation, is J. Jeans’ treatment (1940) of 
the “persistence of velocity.” 

The problem of non-independence of steps has been treated by various 
authors with different limiting conditions, all of which include the fact 
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that there is xo external force or anisotropy of the surroundings. One of 
their methods assumes that the total dependence of the steps of the ran- 
dom walk is due to a dependence on the direction the particle travels. 
Then, by considering the correlations between the directions, the resultant 
modification of the parameters of the Fokker-Planck equation may be 
found. This was first done by R. Furth (1920), who considered one-step 
dependence and applied the results to a study of the movement of Infuso- 
ria. It was later independently redone by G. I. Taylor (1922) and F. Zer- 
nike (1928), but it was not until 1948 that it was shown (Moran, 1948), 
as assumed by earlier workers, that the introduction of this dependence 
does not alter the form of the Fokker-Planck equation, but merely 
changes one of its parameters. This treatment has been generalized by 
C. M. Tchen (1952) to the case where the dependence need not be one- 
step, i.e., where the partial correlation between steps is not zero. 

The case where the total interdependence of the steps of the random 
walk is due to a dependence on the direction the particle travels may also 
be treated by other means. When there is a constant angle between steps, 
the problem may be solved by the use of coordinate transformation 
matrices, as was done by H. Eyring (1932). This has been further ex- 
tended by W. J. Taylor (1947), H. Kuhn (1947), and M. H. Benoit (1947) 
to the case where the azimuthal angle between steps is not arbitrary. The 
problem of the particles having a finite probability of remaining motion- 
less during the time of a step has been solved by G. Klein (1950-51) by the 
use of the calculus of finite differences. For the case where the lengths of 
the steps and the time between the steps can take on only one value and 
the particles are confined to move along a tetrahedral lattice, G. W. King 
(1951) has re-derived the properly modified Fokker-Planck equation di- 
rectly, also by the use of difference equations. 

For the case in which not only the directions but also the speeds of the 
walk are correlated, W. H. Furry and P. H. Pitkanen (1951) have derived 
a solution for a specific application to gaseous diffusion. 

These extensions of the random walk method have been applied suc- 
cessfully to the case of gaseous diffusion (Furry and Pitkanen, 1950; Yang, 
1949) and to the long-chain molecule considerations (Guth and Mark, 
1934; Taylor, 1948). 

In this paper we shall consider the situation in which only the directions 
which the particle travels are correlated, and then only one-step correla- 
tions, i.e., the partial correlation between nonsuccessive steps is zero when 
the effect of the intermediate steps is eliminated. However, we will also 
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assume that there are external forces and anisotropy of the surroundings. 
Thus, although in certain respects we will have a more limited treatment 
than has been given hitherto, in others we will have a more general one. 
The results of this paper will thus be applicable to certain physical situa- 
tions where previous treatments were not suitable. 

An example of a possible application of this paper can be found in the 
genetics problem mentioned previously. In this problem the parameters of 
the Fokker-Planck equation are not constants and hence the earlier meth- 
ods are not applicable. If any of the three causes for systematic change in 
gene frequency (Wright, 1945) is considered to be dependent not only 
upon the present state of the population but also upon its previous state, 
then the formulas of this paper [eq. (43)] may possibly be used. 

Another situation in which the results of this paper may be applied is in 
diffusion with an external force or in an anisotropic medium. Also, in the 
study of long-chain molecules under an external force and with the “lining- 
up” of other long-chain molecules in the medium about the molecule under 
investigation, the formulas derived in this paper will perhaps be useful if 
we assume that to a first approximation all azimuthal angles are equally 
probable (see the second example of distributions). For the motion of or- 
ganisms under an external stimulus, the results of this paper will be ap- 
plied in a future article by the author (1953). 

Definitions and assumptions. Instead of speaking of the probability that 
a particle is at a point, it is conceptually easier to speak of a large number 
of particles moving around, the density of the particles about a point 
being a measure of the required probability. In the usual random walk 
picture one speaks of a particle traveling a certain distance—there being a 
certain distribution of these distances—and then turning. In our picture 
we will speak of a particle traveling in a straight line for a certain length of 
time 7 with an average speed c before turning, where turning means a 
change in direction of the particle’s motion. Pearson (1906) and Yang 
(1949) have called 7 the “time of flight.” Since we will not be interested in 
a particle’s speed fluctuations during its travel between consecutive turns, 
the speed c of the particle will be taken as the average speed of the particle 
between turns. There will, of course, be a distribution of the 7’s and the 
c’s, these distributions playing roles similar to the distributions for speed 
and distance that are commonly used in kinetic theory discussions. For 
ease of discussion we may consider the particles in anthropomorphic 
terms, i.e., the particles may be assumed to “know” what their c and 7 are 
at the start of their “trip” between two successive turns. 

In order to make the idea of a random walk completely explicit—as 
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opposed to diffusion—we will assume that the particles have negligible 
interactions with each other, and thus collisions between the particles in 
question can be neglected. 

The following assumptions will be used throughout the paper: 


1. The particles have negligible interaction with each other. 

2. Each time the particle turns it starts off afresh, with no “memory” 
of its previous ¢ or T. 

J. The amount of time spent in turning is negligible compared to the 
time the particle spends traveling between turns. 

4. During a unit length of time the number of particles in each small 
region, as well as the distributions of c and r, remain approximately the 
same (i.e., the system changes slowly). 

5. All the functions in this paper change sufficiently slowly over the 
average distance the particle travels between turns to allow us to drop all 
higher terms in a series expansion. 

6. All the functions in this paper are continuous and analytic. 


We shall first derive the net “flow” of particles across a given surface by 
a method similar to that first given by J. C. Maxwell (1927) in his deriva- 
tion of the diffusion equation of a gas and given in greater detail by 
M. Loeb (1934). This derivation will be carried out under the assumption 
that both the distance traveled and the probability of a given direction of 
travel vary with the direction. Then it will be shown how the problem of 
persistence and external bias, which will be defined later, may be handled. 
After this is done certain specific distributions will be introduced, the re- © 
sults of which may be compared with results derived by other methods. 
We will also illustrate how these results may be applied to specific situa- 
tions. 

Flow equation. We shall try to find the net number of particles 
Tro, yo, 20 (2, £)dxdy that, per unit of time about the time é, crosses the plane of 
dimensions dxdy located perpendicular to the z-axis at the point (%o, Vo, Z0), 
where by net number we mean the difference between the number of par- 
ticles crossing the plane from the region z < 2 and the number crossing 
the plane from the region z > 20. To find Vx, wo, (z, t)dxdy, the follow- 
ing program will be adopted: a) find the number of particles with an 
average speed between turns between c and c + dc, a time of travel be- 
tween turns between 7 and 7 + dr, traveling in the direction toward the 
point (xo, Yo, 20), and which turned in the region between « and x + dx, y 
and y + dy, and z and z + dz during the unit time about ¢; ) multiply 
this by the probability that the particles will reach or pass the plane 
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dxdy at the point (%o, yo, 20) before making a subsequent turn; c) sum 
(by integration) separately over all the positions, c’s and 7’s, for z < 20, 
and for z > Zo; and d) take the difference between the two sums. [Here- 
after when it is stated that a particle has a given value of a variable (e.g., 
c), it will be understood that this actually means that the particle has a 
value of the variable within a small range of the variable (i.e., in the 
region between c¢ and c + dc).] It will be found convenient in certain 


(XoYoZo 


FIGURE 1 


places to express the point (x, y, z) in terms of spherical coordinates 
(r, 0, 6), where the origin of the spherical coordinates will be taken as the 
point (x0, Yo, Zo). (See Fig. 1.) To simplify the notation the point (x, y, z) 
will be represented by the symbol X and, analogously, the point (ao, yo, Zo) 
by X 0. 

Let the number of particles in the region about the point X at the time ¢ 
be N(X, t)dxdydz. Let P(c, r, 01, ¢1 |X, é) sin 0,d¢,d0,dcdr be the probability 
that a particle has an average speed and time of travel between turns of c 
and r respectively and is traveling in the direction given by the angles 
(01, $1), provided that it is at the point X at time ¢. Fora particle which is 
at the point X to be traveling in the direction toward the point Xo in the 
plane dady it must be traveling in the direction given by the angles 
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(x — 0, r + ¢). The solid angle sin 6:d¢,d6, is in this case given by 
|cos 6 |dxdy/r?. Then the number of particles with speed c and travel time 
7 traveling in the direction toward the point Xo in the plane dxdy and 
which, during a unit time about ¢, turned in the region about the point X 
is given by the product of N(X, #)r* sin 6drdgdé, 


P(c,7,7—0,7+9|X, 1)( [0s | a cdr, 
and the number of times a particle starts a new path (i.e., turns) per unit 
time. Since 7 represents the time that a particle travels between turns, and 
since from assumption 4 we have a quasi-steady state, the number of times 
that a particle turns per unit time is 1/7. Then the required product is 
N (X,t)P(c¢, r, r— 0, 7+¢|X, t) 


T 


x 


|cos 6 | dxdy (1) 


r2 


r2sin @drd¢d6dcdr. 


The probability that a particle which turns at a distance r from the 
point Xo with a given c and 7 and is traveling in the direction toward the 
point Xo will reach this point before turning is, from the definition of c 
and 7, obviously 1 if r < cr and O if r > cr. Hence since cos @ is positive 
for 0 < 6 < 7/2 (ie., for z > 20) and negative for 7/2 < 6 < x (ie., for 
z < 29), therefore 
Tx, (2, t) dxdy = 


fff ff 2S )P(c, 7, See Lean ) (9) 


yg £08 Lees 


r2sin Odrdgodédcdr. 


If the integrand of equation (2) is expanded in a Taylor series about the 
point Xo, keeping only the terms through the first derivative (because of 
assumption 5), we have, to a first approximation, 


Ix, (2; t) ned 


a Ae A A a Ae kee NPC +, eH 0, FPS XD) 


+r sin 0 cos @ IN (X, )P(c 1) 7 8, +9] X51] b 


+rsin 0 sin 6 oN (X, t)P(c, 7, 7— 0, r+ $|X, t)] 


+ 7.cos oS IN (X, t)P(c, tr, 7— 9, r+¢|X, t)] by 
Xcos @sin 06 drd¢dédcdr. 
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In order to further evaluate the above expression, it will be advanta- 
geous to find a more useful expression for P(c, 7, 61, ¢1 |X, ¢). This will be 
done in terms of the probability p(c, 7, 01, ¢: |X, £) sin 0,d¢d0,dcdr that a 
particle that turns at the point X at time ¢ will then “acquire” a speed c, 
travel time 7, and direction (;, $1). Since a quasi-steady state has been 
assumed (assumption 4), the number of particles with a given c, 7, 6:, and 
¢: which turn in a region about the point X in unit time about ¢ will be 
equal to the number of particles which, after turning at the point X in 
unit time about #, will acquire the same c, 7, 6, and ¢;. The former is given 
by the same argument as was used to derive equation (1) and therefore is 

N(X,1)P(c, 7, 01, ¢1| X, 2) 


T 


sin 0,d¢,d6,dcdrdxdydz. (4) 


Since the number of particles that acquire a given c, 7, 6,, and ¢; ata 
point X in unit time about ¢ is equal to the total number of particles 
which turn at the point in unit time about #, multiplied by the probability 
that a particle that turns at the point X at time ¢ will then acquire a given 
c, T, i, and ¢1, we have 


INA Kb EG Es, HG; ak 
Sia abt ttbndrh Xa lee toean 61, b1| X, #) 


co foo] . aA CX Pie , 1, ei 
sp { fPae : es ee Ae 6,d¢,d6,dcdr. 


0 0 T 


Since P(c, 7, 1, gi |X, t) isa probability, the integral over the whole range 
of its variables is equal to 1, and therefore 


PUG, T, 4;, $1| X, t) 


tp (ec, Ty) 61, o1| X, t) (5) 
fo) oo © 2a ¥ 
ifs ule of re rp (c, 7, 0, i] X, #) sin 0:4¢140.d cdr 


0 0 


Let the bias probability, B(@, 1 |X, é) sin 9d¢id1, be the probability 
that a particle will travel in the direction given by the angles (61, ¢), if 
it turns at the point X at time /; let pilc, 7 |1, di, X, t)dcdr be the prob- 
ability that a particle will have speed ¢ and a travel time T, if it turns at 
the point X at time ¢ and travels in the direction given by the angles 
(#1, 1). Therefore, by definition of conditional probabilities (Cramér, 
1946), 

plc, T, 41, g2| X, t) =B (6, $1| X, t) Pile, T | 61, $1; +e t) : (6) 
It will be assumed that the probability of traveling in any direction has 


only a slight dependence on the direction. The bias probability may there- 
fore be expanded in a series of spherical harmonics and to good approxi- 
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mation only the linear terms need be kept. Thus 


1 
B (61, $1|X, i) = =— 11+ B..CX, t) sin 6; cos $4 
Ar (7) 
+ By (X, #) sin 6; sin ¢1+ 8. (X, t) cos 04), 
where 6,(X, ¢) (the coefficient of bias in the x-direction) is the bias the 
particle has to travel in the positive x-direction, and similarly for 8,(X, £) 
and 6,(X, ¢). Since B(, ¢: |X, #) is always positive, therefore 


Bete et Peaks 8) pst) 1) (8) 


In most problems that have been treated using the random walk meth- 
od, the analogue of fi(c, 7 |6:, ¢1, X, ¢) has been considered to be inde- 
pendent of the direction of travel. However, we shall consider the situa- 
tion where p:(c, 7 |@1, ¢1, X, ¢) is slightly dependent on the direction. This 
could be formulated [in an manner analogous to the expansion of 
B(A1, $1 |X, #)] by saying that 
pile, T | 61, 1, As t) = py(c, T) EX'S t) + ,(¢, Ty xX; t) sin 0; cos ont 

+ py(c, 7, X, t)sin 6, sin dit p.(¢, 7, X, t) cos 41 


| pei 6, ty X, 1) | << pole; T, X5t), = 1, 2,3), 


and that 
where 
t= 5 Xe=y, ty SB se (9) 


However, a less stringent formulation for pi(c, 7 |1, $1, X, 4) will suf- 
fice for this paper. That is, 


ff orm ps (6, + | Ob X, 0) ded 

0 0 
df J c*r™ [py (c, r, X, t) + p2(¢, 7, X, t) sin 01 cos $1 (10) 
+ p,(c, 7, X, t)sin 0: sin dit p.(¢, 7, X, t) cos 61] dedr 


(n=0,1,2; m=1, 2) 


and 


iG <P ae 
ee eee kn ids (6 Sy es bh ado 


where a, b, a’, and 0’ satisfy the relation 
[Of arms (o 1| Ox by X, 1) dede 
0 0 


(12) 
a fl fC orm ple, r | 61, 1) xe t)dcdt 


(n=0, 5 m = 1, 2). 
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That is, inequality (11) is true in the region where pi(c, r| 01, $1X, #) is “im- 
portant.” Since the angles (6, ¢1) are arbitrary, therefore substituting 
equation (10) into equation (12) and equating the coefficients of the con- 
stant term, the sin 4, cos ¢; term, the sin 4, sin g; term, and the cos ¢, 
term, we have 


ff or be (65.% XX, DACar 
0 0 
[lf 22 be, (6, Ts we t) dcdr CLS) 


(n= 051,2 + -m = 12s tei Fhe 


Therefore, from equations (10), (11), and (13), we have 


Sf or ba: (6, tT, X,t)dcdr 
0 ~o 


aw a) uv b 
ee BE Ws c"r™"..( 6, 7, X, Ddcds Sefd f ct; ™ 
a’ Ja s a a 


X | Pa (67,X, 0) |dedr< [P [crm py(c, +, X,t)dedr (14) 


~ ff or pole, tT, X,t)dcdr 
0 ~o 
(s=0, 1,2; m=1,2- = 1 253), 


The following notation will be used: 


Sf Sf orm pale, 7, X, t)d cdr = (F(X, 1 
0 0 


(15) 
JL fl 01" be, (6,1, X,)dede = (cr) 5, (X, 1) (i= 1, 2, 3). 
0 0 
Thus by the use of the relations in (14) and (15) we have: 
| Corry tay t) |< (er) (XX, t) . 
(16) 


(n=0, 1,2; me 1 2. te i ool. 


The expression (c*7™)9(X, £) may physically be considered to be an aver- 
age value of the non-directional component of c*r™ of a particle which 
originates from the point X at the time t. Similarly (e"r™),,(X, )(¢ = 1, 2, 
3) may be considered to be an average value of the directional component 
of c"r™ in the direction «; of particles which originate from the point X 
at time ¢. For example, if n = m = 1, then (¢7)o(X, #) may be considered 
to be the average distance the particle moves independently of direction, 


RANDOM WALK WITH PERSISTENCE AND EXTERNAL BIAS 321 


if it starts out from the point X at the time é; while (7 )ailaX yt) (¢. = 915-273) 
would be the additional distance the same particle travels in the direction 
of positive x; over and above its average distance of travel. 

A final point to note is that if f(c, 7, X, #) is an arbitrary function of 6; 
7, X, and ¢ and obeys assumption 6, then 


[Of orm Uo, ck, Nedra ff orm 


Xf (¢, t; X, i) dcdr ($71,125.55). 


Hence inserting the expressions of equations (5), (6), (7), and (10) into 
the right-hand side of equation (3), carrying out the integration over r, 0, 
¢, c, and 7, and using equations (8) and (16) to simplify the results, we 
then have to a first approximation: 

(er) ,(X,#) , Bs(X, t) (er) 0 (X, 2) 


ee 


Dixy (2,1) ={[S “He oe ils SU a) 


pans ‘) 
(17) 


To ate 1) 
' By the use of the familiar argument of considering the difference of net 
flows across the faces of a cube of size dxdydz, we find that 


ON (X, t) 33 — [Sx t) or'x (y, t) 4 0x (2; 24 


— (18) 
Ot Ox Oy 02 


where I'y(x, t) and I'x(y, é) are found in a manner similar to that used for 
finding I'x(z, #). Therefore, substituting (17) into the above equation, we 
have, V being the ‘‘del”’ operator, 


272) 9 (X, sear ene 
gh abd 


ath! 


= 21 N 
al 3 aA 


oe 
(19) 


(cr) = (X, t) eee! t) (er) 0(X, 2) It. 

79 (X, t) 79(X, #) 

In S. Chandrasekhar’s berinolony (1943), < x; > represents the net 
average distance and < x; > represents the average of the square of the 
distance, both in the x; direction, that the particle travels between turns. 
From the description of what (cr)a(X, 1), Bx(X, 4), (cr)o(X, #), and 


(€273) 9(X, #) signify and from equation (16), we see that 


x [ wcx, t) 


(er) 2, (X, t) + Be, (X, t) (cr) (Xs t) = <ai> , 


(Ar X DS <aIS 
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Therefore, we see that (19) is approximately equivalent to the Fokker- 
Planck equation [see Chandrasekhar, Joc. cit., eq. (126)], as would be ex- 
pected. The term “approximately equivalent” is used since we are con- 
sidering a non-zero average time between turns, which is slightly different 
from that of A. Kolmogoroff (1931), who considered the limit of his aver- 
ages for the travel time going to zero. If the assumption given by (16) is 
not used, then the resultant differential equation may also be shown to be 
approximately equivalent to the Fokker-Planck equation. 

Flow equation with persistence of direction and external bias. A coef- 
ficient of bias, 8,,(X, ¢), has been introduced to account for the fact that 
the probability of a particle traveling in a given direction after turning 
may not be the same for all directions. The bias may be divided into two 
parts: an external and an internal bias. The external bias is that which is 
due to forces which lie external to the particle, i.e., variations in the me- 
dium that the particle is traveling in, external force fields which may 
affect the particle, etc. Also, if the particle has a preference of turning 
toward a given direction regardless of the direction it was traveling before 
it turned, that may also be classed under external bias. It will be explicitly 
assumed that the effect of this type of bias is independent of the direction 
in which the particle was traveling before it turned. 

The internal bias is that which has been called persistence of direction, 
1.e., the particle tends to persist in traveling in the same direction which it 
was traveling before it turned (Jeans, 1940). This may also be stated by 
saying that the direction the particle travels after turning is dependent 
upon the direction the particle was traveling before it turned. It will be 
further explicitly assumed that these two types of biases are independent 
of each other, and thus the probability that a particle turns toward a 
given direction is the product of the probabilities that it turns toward this 
direction due to the external and to the internal bias. 

Let E(2, ¢2|X, t) sin O.dg2d0, the external bias probability, be the 
probability that a particle, after turning at the point X at time ¢, will, due 
to the external bias, travel in the direction given by the angles (02, ¢2) and 
let [(62, $2 |01, 1, X, é) sin Oxd¢xd6, the internal bias probability, be the 
probability that a particle, after traveling in the direction given by the 
angles (6;, ¢,) and then turning at the point X at time ¢, will, due to per- 
sistence of direction, move in the direction given by the angles (6, do), 
where the angles (41, ¢i, 62, g2) are as shown in F igure 2. Letting 
A (82, 2 |01, $1, X, t) sin 0,d¢d0» be the net probability that a particle after 
traveling in the direction given by the angles (61, ¢:) and then turning at 
the point X at time ¢ will move in the direction given by the angles 
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(02, 2), then, due to the independence of the probabilities, 
A (62, $2 | 61, $1, x t) 


= E (42, da| X, t) I (42, d2| 41, 1, X, t) 
T 2x ’ 
Jf E (Ox 21, 9 1 (G25 $21 81) 15 X, 0) sin Oadead 0g 
0 0 


where the denominator is put in to make it a probability. Let 
B'(61, $1, 82,2 |X ,¢) sin 6 sin 62d¢d0,d¢2d6, be the probability that aparticle 
which turns at the point X at time ¢ will have turned there after traveling 
in the direction given by the angles (6;, ¢:) and will, after turning, travel 
in the direction given by the angles (62, ¢2). This is obviously equal to the 


(20) 


Z 


~ 
~ 


1 
al 
| 

» | 
! 
q 


= INCOMING PATH 
—-— =OUTGOING PATH 


FIGURE 2 


probability that the particle which turns at the point X at time ¢ reached 
this point after traveling in the direction given by the angles (41, 9) times 
the conditional probability that if the particle did turn at the point x at 
time ¢ after traveling in the direction given by the angles (1, 1), it will 
then travel in the direction given by the angles (, ¢:). Hence if 
W (61,1| X,¢) sin 0:d0:d¢ is the probability that aparticlewhich turns at the 
point X at time ¢ will have reached this point after traveling in the direc- 
tion given by the angles (1, ¢1), then 


B’ (01, 1 02, ¢2| X, #) =W (61, d1| X, t) A (42) $2] 01, dy X, 4). 
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Since obviously 
B(6@ p X = a Bal ps t) sin 6,d¢1d 0; 
? ? t ) ? » Y2 5) ; 
( 2 2 | ) if A; i Ao p 


therefore 
B (62, 2, XxX, t) 


(21) 

Le Qa 

= [" f WO, 411 X, 0 A (8, $2| 815 2, X, 0) sin Odd. 61. 
0 0 


In order to evaluate the above integral, several relations must first be 
established. We must first find a more “useful” expression for W(h, 
1 |X, t). If w(1, d1, X, £) sin 0,d¢,d6,dadydz is the number of particles that 
turn, per unit time about the time #, in the region about the point X after 
entering this region from the direction given by the angles (6;, ¢;), then it 
is obvious that 
W (61, $1|X, #) = ATR 88 


7 2a . pl, 
if: w (41, $1; XxX; t) sin 6,d¢,d A, ( ) 
0 0 


For a particle with speed c, travel time 7, and traveling in the direction 
given by the angles (6,, ¢:) to turn in the region about the point X, it must 
previously have turned in the region about a point of distance cr and 
direction (7 — 6, 7 + ¢1) from the point X. Thus if we take the origin of 
the spherical coordinate system as the point X (instead of the point Xo as 
used previously), then the number of particles with a given c, 7, 6:, and 
¢: which enter the unit region about the point X per unit time will be equal 
to the number of particles with the same c, T, 6, and ¢; which leave the unit 
area about the point (cr, 7 — 61, + ¢) per unit time. Then, by the use 
of reasoning similar to that used in the derivation of equation (1), we see 
that the number of particles with a given c, T, 0, and ¢; which leave, per 
unit time, the area about the point (cr, r — Hh, r+ gy) is 

N(¢r, t— 01, r+ 4, t)P(c, 7, 6, $1| ¢r, r— 61, r+ 44, 2) 


T 


Xdcdr sin 0,:d¢,d 0, . 


Thus the total number of particles which enter with a direction given by 
the angles (6, ¢:) and which turn, per unit time, in the unit area about the 
point X’, is obviously the integral over all c’s and 7’s of the above expres- 
sion. 


Therefore 


w (1, d1, X, t) =f Se 


N (er, (eer: 61, t+41, t)P (c, ns 61, $1 


(23) 
a a EAE al Or oO Kee ee 
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We will now explicitly assume that the change of the flow of particles 
across any unit plane in any of the x, y, z directions at any point is much 
less than the number of particles which leave the unit region about the 
point per unit time, i.e., 

Oly (x;,¢t) | N(X, ?#) eee: 23 


=> p) 
bee FD i eek (4s) 


where the right-hand side is found by integrating (4) after inserting equa- 
tions (5), (6), (7), and (10) into it and simplifying by use of (8) and (16). 
From the above and (18), we see that the number of particles which enter 
and leave a region per unit time is approximately the same (which is in 
agreement with assumption 4). If we expand the integrand of (23) in a 
Taylor series about the point X, keeping terms up through the second 
derivative [due to assumption 5 and the fact that we have effectively 
integrated over the variable rv of expression (2)]; insert equations (5), (6), 
(7), and (10); integrate; simplify by use of (8), (16), and (24); and sub- 
stitute the results into (22), we have 


Wee (er) o(X, t) 
aA Caer Toh Xad) 

N (X,t) 

rg t) 


W (0, 6|X,) =7-4 1+] 6. (X, 0) — 


F oy To (X, t) 25 
Xsin 6; cos $1 + | By (X, t) — V(x, 1) (25) 


T9(X, t) 


2{ W(X, 2) Deane 2 
Xsin 6; sin $1 + | B,(X, 0) eR TTS CET RT cos 61 


To (X, t) 


As a check of the above method, if the simplification introduced by (8), 
(16), and (24) is not used in the evaluation of w(H, di, X, #) [except in the 
case of having V(X, t)/7o(X, #) instead of the actual value, which is done 
solely for ease in writing and introduces no error for this discussion], we 


we find that 
7 2a ; N(X, t) ON (X, t) / 
i J w (64, $1, ay ft) sin Bb UA IS ae t) zi at & (26) 


Since the number of particles that enters the unit region about the point 
X per unit time, minus the number that leaves [see eq. (24)] is equal to the 
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net change of the number per unit time, which is V(X, ¢)/dt, we see that 
equation (26) is correct. 

Just as with B(62, ¢2 |X, f), it will be assumed that the external bias is 
“weak,” so that E(02, ¢2|X, 4) may be expanded in a series of spherical 
harmonics and only the linear terms need be kept to a first approximation. 
Thus 


E (Ox ¢2| X, 2) we glite(X, i aine oucomtie 
+e, (X, t) sin 02 sin d2 +e, (X, t) cos O2}, 


(27) 


where e,,(X, ¢) (2 = 1, 2, 3) is the coefficient of external bias in the x,- 
direction and may be physically interpreted in a manner analogous to that 
of B,,(X, #). Since E(62, d2 |X, £) is always positive, therefore 


2 (X, 1) +é2(X, t) +2(X,t) S1. (28) 


It will be further assumed that the internal bias is dependent solely 
upon the angle between the incoming and outgoing path of the particle, 
i.e., if a is the angle between the incoming and outgoing paths, then 


I (62, $2| 041, $1, xX; t) = Teh st, By; (29) 


where a is always taken to be between 0 and 7. The foregoing implies, in 
anthropomorphic terms, that the particle has no right-hand, left-hand, 
up, or down internal bias. We shall define the coefficient of persistence of 
direction, W(X, #), as 


¥(X, 1) =2n f'cosal’(a|X,1) sinada, (30) 
0 


shall assume that 
W(X, t) en, (X,t) «1 (¢= 1, 25 3); (31) 


and that the quantity [1 — W(X, #)] is not too small. From the definition of 
W(X, ¢), we see that the more the particle tends to continue traveling in the 
same direction it was traveling before turning, the larger y(X, #) will be. 
If the particle, on the basis of its internal bias, has no preference for any 
direction [i.e., 7’(a |X, ?) = 1/47)], then it can be seen that ¥(X, £) will be 0. 
Also, it can easily be shown that 
aa Lace Dalat ice Hal 
for all Z’(a |X, 2). 
In order to integrate 


wT Qa 
J ys E (4, $2|X, t) I (02, d2| 61, fi; AG t) sin O2ddod O» , (32) 


which appears in the denominator of the right-hand expression of (21),a 
simple coordinate rotation will be useful. This rotation may be done as 
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follows: First rotate the x, y, z coordinate system about the z-axis pi 
radians and then rotate the resultant coordinate system about its y-axis 
9, radians. This will then yield a new coordinate system whose z-axis lies 
along the direction given by the angles (6,, ¢:). Let the coordinates of a 
point in the initial coordinate system be given by (r, 62, 2) and coordi- 
nates of the same point in the rotated coordinate system be given by 
(r, a, w), where the a is the same as the a of (29) and where we may drop 
the coordinate r since it is the same for both coordinate systems. Thus 
taking the matrices of the two rotations and multiplying them together, 
we have 


Sin @ COS w2 
sin a sin We 


cos a 


cos 6,0 —sin 6; cos ¢; Sin ¢; O Sin 45 COS ds 
-(° 1 0 ( sins cro, 0 | [sn sings} (3.3) 


sin 6; 0 cos 6; 0 0) 1 COS 45 


— sin 1 COS $1 0 sin A sin 2 
sin 6; coS ¢; sin 4; sin 1 cos 44 COS A, 


( cos 6; coS ¢; cos 6; Sing; —sin as (sn 85 cos “| 


If we expand the above equation and take the Jacobian, we find, after 


laborious computation, 
O(a,w2) _ sin A 
8 (82, $2) sin a 


and hence, as expected, 
sin 6.d¢20.=Ssin adw.da . 
If both sides of (33) are multiplied by the inverse of the rotation matrix, 
we have 
sin 62 CoS $2 = cos 6; cos ¢; Sin a COS we — Sin ¢; Sin a Sin we | 
+ sin 0; cos ¢; coS a , 
sin 62 sin ¢2 = coS 6; sin ¢; Sin a COS w2+ COS Pj SiN a SiN we (34) 
+sin 6; Sin ¢; cosa , 
cos 6. = —sin 6; sin a CoS w2+ cos 6; cosa. 
Then inserting (27) and (29) into (32), and using (34), we find the value of 
the integral of (32) to be: 


ae, 1+y(X, De (X, #) sin 0: cos di + ¥ (X, t) ey (X, #) sin 41 sin 41 


as +y (X, t)e,(X, t) cos 6}. 
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If we then insert the results of equations (25), (27), and (29) into (21), 
integrate with the use of the coordinate rotation given by (33) and one 
like it except for interchanging the subscripts i and ; [i.e., for integrating 
over the angles (6, ¢1)], simplify by the use of equation (31) and assump- 
tion 5, and insert (7), we have, after simple rearrangements, 


Be, (X,#) =4 [14+-4(X, De, (X, 0 


a (er) 0(X, 2) aes 
CT) 90 ’ 
ow OH) soll Ea ed a 
1 (X,0) WOOD Cnt dBea 
7 1e.¢ t) 


Then inserting this value for 6,,(X, £) into (19), we have, after elaborate 
rearrangements, 


at ween) 


aa 272 
See: =) v( wis By 


“SSE Beet 


(c?r?) 9] \ (er) o 
me ( c CT) =| To | - : 
where it is understood that all terms are functions of position and time. 


If (¢7).,(X, t) and e,,(X, t) are considered to be the x,th component of 
the vectors (cT)(X, t) and €(X, #) respectively, then (36) may be written: 


ae | 


aN 11 i+u( ane ) v(v (ore) 


i) 


(37) 
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If / represents the distance the particle travels between turns, then 
1 = cr, as mentioned previously. Thus instead of using the variables c and 
7, the variables 7 and / or the variables c and / could have been used 
equally well. It may be easily verified that for the former situation the 
expressions (¢7)o(X, ¢), (€272)o(X, #), and (é7),,(X, #) in the above equa- 
tions would be replaced by 7,(X, #), 2(X, #), and ],,(X, 4) respectively, 
where the significance of the terms is self-evident. For the latter case, the 
variables being c and /, the expression (€7)o(X, t), (@272)o(X, #), (€7)z,(X, 8), 
and 79(7X, ¢) in the above equations would be replaced by Jo(X, #), 2(X, 2), 
Iz,(X, t), and (J7¢)o(X, #) respectively. 

Discussion of equations. It should be pointed out that the external force, 
if any, may not only affect the external bias coefficient [e,,(X, ¢)], but may 
also affect the distance the particle travels between turns [i.e., affect 
(cr).,(X, t)]. The methods outlined above do not seem applicable to this 
author if the direction of the particle depends not only upon its previous 
direction but also upon its direction of travel two or more turns previously 
(as occurs, for example, in considering lengths of long-chain polymers). 
(See Taylor, 1948 for a discussion of this point.) 

If 

¥ Xt) a (X= ber) (X94) = 0 tn, 2,3) 4 B38) 


then, from the discussion in the introduction, (36) may be considered an 
“exact” expression of Fick’s second law as applied to a dilute solution. 
Thus the diffusion coefficient D is given by 


_ 1 (c#r4).(X, t) (39) 
ae SEE TE ZO ales 
In this case, (36) reduces to 
eee ee v?[DN(X, 2)] (40) 


as opposed to the usual formulation of Fick’s second law (Jost, 1952): 


aN (X, t) 


+ =V-DV[N(X,#)]. (41) 


If D = D(X, #), then in the steady state in a closed container the usual 
formulation of Fick’s law would predict that since 


N (X, t) =constant , 


which does not appear physically plausible, since for D not to be a con- 
stant there must be variations within the medium and hence the same 
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number of particles at each point would not be expected in the steady 
state. However, modification of Fick’s law, given here, would yield 
constant 

D(X, t) 

for the steady state. It should be noted again that (40) is the “modified” 
Fokker-Planck equation. This equation has been previously derived by 
S. Chapman (1928) for the case of particles undergoing Brownian motion 
in a non-uniform field. It may also be noted that since 


V(X, t) D(X, 1) ] = D(X, NVIN(X,)1+N(X,)V(D(X, 4], 


N(X,t) = 


and if (40) is valid, we may consider the first term on the right of the 
above equation to represent the usual flow term of Fick’s first law, while 
the second term represents the results of an “effective” force due to the 
nonhomogeneity of the medium. 

Equation (40) has been tested experimentally for the case of the Soret 
effect, and has been found to be invalid (Hartley, 1931a). However, 
Hartley attributes this to the fact that the effect of the solvent’s motion 
was ignored in the derivation of (40), and hence it does not apply in this 
case. Hartley (1931b) also attempted to test the validity of (40) by allow- 
ing a second solute to diffuse steadily through a solution, thus making a 
steady non-uniform medium in which the diffusion coefficient for the first 
solute is not uniform. However, the diffusion effect was masked by 
the difference of solubility of the first solute in the different parts of the 
solution containing the second solute, and thus no conclusion as to the 
validity of (40) could be drawn. The latter factor of change in solubility, 
if the diffusion coefficient changes, appears to be generally true for most 
experiments on diffusion in a liquid. This problem might be treated by 
considering that there is an external bias at each boundary between re- 
gions with differing solubilities, such that a particle would “prefer” the 
region in which it is more soluble, in which case the equations in this paper 
might be used. However, this method does not seem very practical. Thus 
at the present time there are no experimental data available for making a 
choice between the usual formulation of Fick’s law or our modification. It 
should, of course, be emphasized that (40) was derived for particles under- 
going a random walk, which is certainly not the same as particles diffusing 
in a liquid. 

If D is a function of V, then in the steady state we would have 

N (X, t) = constant 


for the closed container for both formulations of Fick’s law. However, for 
D to be a function of V, assumption / would have to be violated. If D is 
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only slightly dependent on N, then it might be expected that the diffusion 
equation would be approximately the same as for the case where D is a 
function of X and ¢ only, and therefore (40) may be true for this case, also 

For the case where the particle is restricted to movement in a plane 
(random walk in two-dimensions), the method of treatment, assumptions 
made, and results (except for the constant factors) are completely analo- 
gous. That is, if Y and all other functions are used in place of (x, y) and if 


Qa 
WY D= f cosa I’(a| Y,¢t)da, 
0 


where J’(a |Y, ¢) is defined in a manner similar to their three-dimensional 
counterparts, we find that 


aN_1v a ? i+¥( S- ) 3 (x Pe) 


Ot) 24 dx, )2 1—y ax; 


1 


To 


a [i S4+((1+¥le, (42) 


wv (er)o 9 [(c?r*)07\ (er) 0 
IE rane) Gehan Go| To |y y 


where it is understood that all terms are functions of position and time. 

For the case where the particle is restricted to movement on a line 
(random walk in one-dimension), the method of treatment is again com- 
pletely analogous to that of the three-dimensional case. However, the 
assumption of a ‘“‘weak’’ bias need not be made. This is because, if the 
functions are defined in a way analogous to those of the three-dimensional 
case, and if we let the (+) sign be for the particle traveling in the direction 
toward increasing « and the (—) sign for decreasing x, then B(+ |x, t) and 
E(+ |x, #) may be broken up into: 


B(+ |x, #) et 


(43) 


exactly. - 
The coefficient of persistence (x, ¢) is defined such that the probability 


of a particle maintaining the same direction after “turning” at the point 
x at time ¢ due to persistence is 3[1 + (a, #)], while the probability of the 
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particle traveling in the opposite direction after “turning” is 3[1 


— (x, t)]. Then, if the assumption that (x, f)e-(x, t) « 1 is not made, 
we find that 


2 Caren 
aN @]1 itv (1d [? (c272) g 1] ve) 
“Ot Ox \2 1—y[1—é] —yé 


where it is understood that all terms are functions of position and time. 
If it is now assumed that 


then (44) reduces to: 


on te! VC cee, =) 2 (wie) 


To 


- [S24 (tye 40) 


T=¥ Cor bel Cone] IY}: 


which is completely analogous to those before. [It should be reemphasized 
that the assumption of (45) is (x, t) «1, not W(x, f)e,(x, t) K 1, as one 
would expect it to be in analogy with the two- and three-dimensional 
cases. ] 

Examples of distributions. In order to illustrate the rationale of some of 
the assumptions made, two examples will be given here. 

1. Diffusion. In the considerations of diffusion of solute particles in a 
solution or of Brownian motion, the particles are conceived to have some 


—— rf. | 
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speed distribution (usually Maxwellian) which is independent of the direc- 
tion which the particle travels, and which may be denoted by C(e |X, t)dc. 
The distance the particle travels (/) between turns is usually-considered tc 
be independent of direction, and, for the case where the probability of a 
particle suffering a turn (collision) is the same for each point, the prob- 
ability of the particle traveling the distance / may be given by 
(1/)exp(—//X)di, where \ is the mean free path. However, a more gen- 
eral case may be considered using the above formulation, in which the 
mean-free path \ may have a slight dependence on the direction the par- 
ticle travels. For this case, recalling the discussion leading to (7), we have 
A= Ap (X, t) +A, (X, 2) sin 0 cosd ng 
+r, (X, ¢) sin 6sin@d+,(X,¢t) cos 6. ide 
The Ao(X, t) [which is the same as Jo(X, ¢)] in the last expression may be 
conceived to be the “average”’ distance the particle will travel independ- 
ent of direction, while the },,(X, #) (¢ = 1, 2, 3) [which is the same as 
Ie,(X, t)] are the “extra” average distances the particle travels in the x, 
direction. As is true in most physical situations, it will be assumed here 


that 
| Az, (X, t) | Ko (X, 4) (i=1,2,3). (48) 


Thus | 
pi (c,1| 0,¢, X,t)dcdl=C(c|X, t) Sed cdl (49) 


where the is given by (47). If the exponential is then expanded out in a 
power series, equation (49) becomes to a first approximation 


c — 
pi(c, L| 6, >, X, t) =C(clX, Nyy He (Po, O14 4 | 
1 A2(X, 2) .. Mv(Xs t) oa 50) 
bees t) — [2% t) ssi : Lis Ae eI @ t) sin sin @ ( 
Az (X, b) 
+x OL} 


Hence the analogue of (10), (11), and (12), for the case where the variables 
are c and J, is valid for this distribution. We may note that for the distribu- 
tion of (50), 72(X, t) = 2i2(X, t). Thus if 


[ec (eX, Nac A(X, 
en 


then, substituting averages from (50) into the analogue of (36), we find, 
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where it is understood that all terms are functions of position and time, 
oN 1 1 ( 2)-1> {ds 
apr at yl Mere 32 Bn Xo t/t 


+( 1 4W) et yey ay ial Ve|yt- 


(51) 


If the distribution given in (50) is accepted as being correct, we find that 
the diffusion constant of (39) is 


This is analogous to, but slightly different from, that usually found for the 
diffusion constant. If 1/(1/c)(X, #) is equal to Ke(X, i), where K is a 
constant, then the diffusion constant of (52) and that which is commonly 
used would differ only in the value of the numerical constant, which is 
usually of the order of unity. For example, for the Maxwellian distribu- 
tion, K = 7/4. 

2. Long-chain polymers. If long-chain molecules are treated by means 
of the random walk method, as discussed in the introduction, the particles 
may be considered to have only a single time of travel r*(X, ¢). Thus the dis- 
tribution of 7’s may be given by the Dirac delta function, 6[7 — r*(X, #)]dr. 
(This may be used in a rigorous fashion if we consider the Dirac 
delta function as a limit of an error function, for example, but we shall not 
go into that here.) If the particle travels a distance which is always equal 
to the mean-free path d (i.e., \ is the distance between neighboring atoms), 
then the distribution of l’s may be given by the function 6(/ — )di, 
where ) is given by (47). Thus 


pi (l, 710, $, X, t) dldr = 6 [r —7*(X, #)] 5 [1—Al dldr, 


which cannot be expanded. However, it is easily seen that 


i J mrm py (1, 7 | 0, @, X, t) dldr 
0 


it ad ea eT rey ef Az (X, t) 
; Ny CX, Di. : (XD 
X sin @ cos @ tm sin @ sing +m 2 cos gaia. 


Although this equation is slightly different from the analogue of (10), for 
the variables / and r, it will, due to the approximation of the analogue of 
(16), which also holds for this distribution due to (48), yield the same 
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results as before [i.e., the analogue of (36)]. Thus using this distribution 
and substituting into the analogue of (36), we find that, where it is under- 
stood that all terms are functions of position and time, 


“inact eestabd Gules) me? le 


(52) 
ya Xo 
g +((1+¥le, opie Ss] at. 
he; (X, #) =e,,(X,/) =0, 
vy (X, t) =W= constant, 
HOE t) = Ao = constant , 
7*(X,t) =7*= constant, 
then (52) reduces to 
aN(X,t) _1+¥ % 
Seas ~Toyp 672 VINX, 1, (53) 


which agrees with the equations derived by other methods (see introduc- 
tion for references). 

For the case of a long-chain polymer composed of carbon atoms, the 
angle a between the incoming and outgoing directions is (Fig. 2) 
(a — 109°28’). Thus 


I’ (a| X, t) sin ada = 6 (a — 109° 28’) da , 


and therefore 
y (X, t) = cos (r— 109° 28’) =2. 


Hence equation (52) reduces to 


A3Y [vs] - sLas ES 


(54) 
+ (Feit +5 [dol )*]rt 


If we consider the case where r*(X, #), Ao(X, #), €2,(X, #), and Az,(X, 2) 
are constant, then using the vector notation of (37), (54) reduces to: 


tant A fe NmretheeXe |e oy ee 5) 
< ptvin(x, 4) —5[+3e75 VIN(X,#)]. ( 


Since ¢ is equal to mr*, where is the number of bonds in the molecule, 
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then we can easily find from (55) that the average length of the long-chain 
molecule is 


16 1/2 
blatgeroin+(5° uD) : (56) 


If we consider that the external force is along the x-axis and a non- 
isotropy of the medium is introduced due to all of the molecules “‘lining- 


up”’ in this direction, then i 
In| =, 


je| =e,. 


Thus the average length of the long-chain molecule will be 


1/2 
FO+$eh) + (Gl) dor, (57) 


a result which may be applied to discussions of entropy changes in long- 
chain molecules. 
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In a previous paper, in which a possible mathematical approach to history was outlined, it 
was shown that urbanization plays an important part in the propagation of new ideas. The 
rate of such propagation influences the rate of historical developments. The present paper 
deals in more detail with possible mechanisms of formation of earliest cities. Equations are 
derived which give the limiting size of such cities and their rate of growth. Of particular im- 
portance for the spread of new ideas is the spread of information. The latter largely depends on 
the fraction of individuals who travel between city and country. Expressions for this quantity 
are derived. An approach is outlined to the mathematical study of the earliest social classes, 
which may have been formed as a result of military, religious, or economic stratifications. 


In a previous paper (Rashevsky, 1953), in which we outlined a possible 
mathematical approach to history, we pointed out the importance of ur- 
banization as a factor determining the rate of historical developments. 
Society changes and develops as new ideas, originating with a very small 
minority, gradually spread through it. The society can adopt new ideas 
and new behavior patterns only to the extent to which it is informed about 
those new ideas. In sparsely populated rural communities the spread of 
information is very slow. In cities it is much faster; hence the importance 
of urbanization as a determinant of the rate of progress. 

In this paper we shall further develop some simple mathematical ex- 
pressions pertaining to sizes and number of early cities. We shall also dis- 
cuss some other problems of history, which suggest possibilities for a 
mathematical approach. 

One thing must be very strongly emphasized at this point. We are at- 
tempting here the very first steps of the application of mathematical 
methods to history. Therefore, our immediate aim is not so much to de- 
velop some generally applicable theory as to develop the methodology of 
the new approach. We shall do this by discussing some more or less 
plausible situations in regard to the mechanism of formation and growth 
of early cities, or the mechanism of formation of social and economic 
classes in early society. However those situations need not necessarily 
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have actually occurred. A critic, who might remark that things actually 
did not happen in such a way according to present day evidence, would 
miss the point. The situations are sometimes chosen more for their con- 
venience of illustrating the method than for any other reason. A mathe- 
matical scientist will see that even if the situations discussed here as ex- 
amples did not actually occur, the same method of approach can be ap- 
plied to any situation which the historian or archaeologist may find much 
more likely. We do not yet propose to develop a systematic theory of 
specific historical phenomena here; all we want to do is to show how his- 
torical phenomena can be discussed mathematically. If, by chance, some 
of the situations discussed here turn out later to be incorporated in a con- 
sistent theory, so much the better. But this is not presupposed here. All 
the above must be kept in mind by the reader. 

1. Formation of villages. We shall attempt here to give some mathe- 
matical models of the gradual formation of cities in an originally peasant 
village culture. Closely related to that problem is the formation of early 
economic and social classes. 

Paleolithic man must have been a rather rare animal, and his popula- 
tion density must have been small. Our very tentative estimate (Rashev- 
sky, 1953) of 1 individual per 10 km.’ is not in discord with A. L. Kroeber’s 
(1939) findings for pre-Columbian American Indians. The ‘‘food-gather- 
ing” man must have lived in small groups, perhaps of one or a few families. 
During the nomadic stage the splitting of the group into independent units 
may have been easier than at later stages, when mansettled down to primi- 
tive food production. With small population densities, land must have 
been plentiful, except for arid regions of desert-like nature. A family or 
group of families, once settled down, tend to remain together, both from 
the sense of gregariousness and for the practical purpose of defense against 
enemies of all kinds, whether other men or animals. How large would such 
a group grow? 

The smallest crop yields now in undeveloped countries, taking wheat 
as an example, are about 300 kilograms per year per hectare (French 
Morocco), and this is an unusually low figure (Statistisches Reichsamt, 
1937, p. 42*, ff.). The world average is 970 (without the U.S.S.R.) and 
Western Europe has yields as high as 2500 to 3000. But even 300 kilo- 
grams a year are approximately sufficient in caloric value to maintain one 
individual. Even with the primitive methods of the earliest peasant village 
cultures, yields in the fertile valleys of the Nile, Tigris-Euphrates, Indus, 
and Huanho must have been much higher. However, using the lowest 
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figure, we find that 1 km. of arable land will give subsistence to approxi- 
mately 100 persons. Considering that not all available land is arable, that 
lack of knowledge of crop rotations made it necessary for one person to 
shift from patch to patch, and that grass land was necessary to provide 
for a minimum of domestic animals such as horses or cows, actually much 
more than 1 km.? of land will be needed to support 100 persons. With suf- 
ficient safety we may multiply the amount of land by 10, giving us about 
10 km.’, for a group of 100 persons. This area corresponds approximately 
to a circle of radius of 2 km. Thus a village with 100 inhabitants would 
have to cultivate land within a radius of about 2 km. This may already be 
close to the limiting distance for those who have their lots near the periph- 
ery. Putting that limit at 5 km., we find an area of 75 km.”, which is enough 
to support a village of about 750 inhabitants. This figure gives us the ap- 
proximate upper limit of a village. Beyond that size either single families 
or, for reasons of gregariousness and especially safety, larger groupsof fam- 
ilies will have to migrate beyond the bounds of the village proper and form 
new villages. It is interesting to note that early American Indian villages 
ranged in their populations from one to seven hundred individuals (Swan- 
ton, 1946). The ancient village of Jarmo had about 200 inhabitants. 

Even if the density of population in this early neolithic stage has in- 
creased, it still must have been small. If the average population density 
was 1 ind. km.~?, then, taking 500 individuals as the average size of a 
village, we find one village per 500 km.? This gives us an average distance 
between the villages of the order of 30-40 km., a distance found to hold for 
earliest nome capitals of Egypt (Petrie, 1923, p. 4). 

2. Formation of classes. What can cause the appearance of social and/or 
economic classes in such an agricultural community? 

As we have shown previously (Rashevsky, 1948, 1951) socioeconomic 
classes may be formed on the basis of individuals tending to associate only 
with other individuals who are sufficiently “‘similar’”’ to them in some re- 
spects. Hence a dissimilarity of some kind must be present to cause a 
stratification of a society. Such dissimilarities may be of great variety 
Perhaps the most important are those associated with military leadership, 
religious leadership, and economic status. re 

Inequalities of status are inherent in any military organization, no 
matter how primitive. The primitive chieftain or the modern general is 
distinctly different from the rank and file soldier. Thus, from the earliest 
times, a military class was formed. The core of that class is made up of the 
military leaders and chieftains. The ratio of the number of the leaders to 
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the number of the rank and file soldiers seems to be fairly constant under 
various conditions and is of the order of magnitude of 1%. Thus, even in 
the extreme case of a society in which every male is a soldier, the size of the 
ruling military class would be of the order of 1% of the whole population. 
In less militarized societies it will be even less. Thus, assuming the above 
computed sizes for villages, we expect one or at most 2 or 3 military chiefs 
per village. The importance and influence of these chieftains fluctuates be- 
tween the times of war and peace, being the strongest in wartime when the 
exigencies of preservation of the tribe may cause everybody to obey the 
chieftain unreservedly. In peacetime the influence may be diminished but, 
with an eye to a possible need for them in a future war, some special sup- 
port will be given the chieftain even in peacetime. The support which the 
chieftain receives and the influence he has can best be measured in terms 
of the economic advantages which he has over the average man. On the 
average, each individual will be willing to support a military chieftain to 
the extent to which he needs him by some form of taxation. In wartime 
the power of the chieftain may permit him to seize all of anyone’s prop- 
erty. In peacetime such a power may be almost reduced to zero. Hence if 
w denotes the average fraction of time which a tribe spends at war, then in 
the first very crude approximation the fraction of the income which an 
individual will give the military class will be equal to w. According to data 
collected by P. Sorokin (1937) and Quincy Wright (1942), the average w 
for all countries does not seem to have varied too much during the known 
span of history, and seems to be of the order of 0.1. It varies, or has varied, 
from culture to culture. In principle, w could be ascertained for different 
cultures, countries, or tribes by historical research. The relative wealth of 
the military class could possibly also be ascertained for various past socie- 
ties. This relative wealth should be equal to w. If a class of the size of 1% 
of the whole population possesses a fraction w of the whole wealth, then 
the average per capita wealth in the military class is 100w times that of 
the average per capita wealth of the population. Whether this is so or not 
could be ascertained in principle. Here is a definite expected quantitative 
relation, challenging the historian to prove or disprove it. 

Somewhat similar relations hold for another important class, that of the 
priesthood, which was very strong in many ancient civilizations. The per- 
centage of priests is again probably of the order of magnitude of 1%. 
Their general and, particularly, economic influence is, however, deter- 
mined by different factors. The priest is the spokesman for the ordinary 
man before various gods, who were more feared than respected by the 
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ancient man (Smith, 1952). In case of a calamity, such as individual or 
mass illness, floods, famines, etc., the gods had to be placated by offerings 
and donations, which were received by the priests as the gods’ official tax 
collectors and treasurers. In times of quiet and prosperity the gods had to 
be propitiated to avert their ire. The priests again were the winners. Thus 
we may expect that the economic and general power of priesthood would 
be directly proportional to the incidence of various calamities. It is gen- 
erally accepted that while the influence of the priests was very strong in 
ancient India, Egypt, and Mesopotamia, it was much weaker in Greece 
and China, but no quantitative figures are given. John A. Wilson (1951, 
pp. 270-71) has estimated that in ancient Egypt the priesthood owned 
or controlled from 10 to 12% of the total wealth and manpower. If such 
estimates can be made for Egypt, it should be possible to make them, 
eventually at least, for other countries; thus giving a quantitative form for 
the above-mentioned generally accepted view. On the other hand, the inci- 
dence of calamities could be obtained either by direct historical statistics, 
or by climatological and geographic considerations. 

A list of historical world famines, floods, frosts, and other calamities has 
been published by C. Walford (1878-79). A glance at the figures shows 
that more than half of all famines listed happened in England, an ob- 
viously absurd result caused by the author’s indiscriminate pooling of the 
readily available data for England with those for other, less well-known 
parts of the world. Thus, for comparative studies this source is of no use. 
Similar compilations, free from the above-mentioned defect, could be 
done, though at the expense of considerable work and time. But such a 
task would be amply justified. 

Preliminary climatological considerations yield the following rather in- 
teresting and possibly significant result. On page 360 of J. T. Trawarta’s 
(1943) book we find a map showing the rainfall variabilities for different 
regions of the world. This rainfall variability is likely to reflect the inci- 
dence of drought and, therefore, famines. Whereas for Northern India we 
find the figure of 20-25%, for Eastern China, where the original Chinese 
civilization was born, we find only 15-20%, and the same figure for Medi- 
terranean Europe. In Mesopotamia the famines were caused not by the 
lack of rainfall, which is not a factor there, but by insufficient flooding of 
the Tigris and Euphrates. But those floods are largely determined by rain- 
falls in the Caucasus, and here the variability is 20-25%. Again in Egypt 
the Nile flood is the important factor. If we assume that this is controlled 
by rainfalls in Central Africa, we come to the figure of 15-20%, the same 
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as for Italy and China, whereas we should expect it to be closer to the 
figure for India and Mesopotamia. There is a discrepancy here. It may 
well be that in the earliest predynastic times, when the original priest 
class in Egypt was formed, the climatic situation was such as to produce a 
greater variability of the Nile floods. 

The discrepancy may, however, have a different explanation. Not only 
is the-frequency of a calamity important, but also its magnitude. The 
effect of less frequent but more disastrous calamities may be greater than 
that of more frequent but less disastrous ones. The effect of a failure of the 
Nile to rise must certainly have been a disaster of most outstanding pro- 
portions. Another important factor which determines the social and eco- 
nomic importance of priests is their ability, through some knowledge 
which is kept secret, to predict the calamity. Such a possibility may have 
been likely both in Egypt and Mesopotamia. 

Another important factor is the different tendencies of different indi- 
viduals to “take chances.” Fundamentally we deal here with the problem 
of risks, though the risks are not calculated, but rather determined by a 
guess or feeling. H. D. Landahl (1951) has developed a theory of neural 
mechanisms which may be responsible for different attitudes in taking 
risks. He finds that individuals can be divided into six different types, 
according to their attitude to taking chances. From further development 
of such studies, it should be possible to derive the relative proportions of 
individuals of different types. These proportions may depend both on bio- 
logical and social parameters. Then from known frequencies of incidence 
of calamities of different intensities, and from the calculated proportions, 
it should be possible to calculate the extent of economic support which a 
given society will give to its priestly class as a sort of “insurance” against 
the ire of the gods. Mutatis mutandis this applies to the military class, 
where we shall have to consider the incidence rates of wars of different 
magnitudes. In this connection the studies of Lewis Richardson (1948) 
may be of importance. Here at any rate is a subject worth a quantitative 
study, for if the quantitative relations suggested here are confirmed, they 
will throw light on the question of why the priest class was stronger in 
some ancient countries than in others. 

The influence of the priestly or the military class may be measured by 
criteria other than their economic status. It is important to emphasize 
this, so as to avoid a possible impression that the suggested possible theory 
is essentially economic in nature. In principle we may measure the impor- 
tance of the military or priestly class by the relative amount of legislation 
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that gives it any kind of privileges, economic or otherwise. Or we may 
measure that influence by the percentage of cases in which priests were 
consulted on various general governmental problems. We could then cor- 
relate those entirely non-economic figures with the frequency of war or 
incidence of natural calamities. The reason why we chose the economic 
status as an illustration is that it seems to be be relatively easier to meas- 
ure; but for historical research the non-economic criteria, put in a quanti- 
tative form, are likely to be of greater interest. 

We have discussed above how a theory could be developed which will 
relate quantitatively social inequalities to economic inequalities. The two 
are usually closely correlated, and it may not be wise to consider one as the 
cause of the other. Both the economic and social status of slaves may have 
originated as the result of formation of the military class, since slaves have 
been recruited mostly from conquered peoples. There may also have been, 
however, purely economic factors, which resulted in the formation of dif- 
ferent degrees of subjugation or bondage of some individuals by others. 
Like in any mathematical approach in the first stage each factor must be 
studied in abstracto separately. In reality probably all factors play a part, 
and the important thing to remember is that the net result of operation of 
several factors is not the same as the sum of the results calculated for each 
factor separately. This is because the interrelations between the different 
factors are most likely non-linear, and the factors therefore are non- 
additive. 

With this idea in mind, we shall discuss here a theoretically possible 
economic factor. We must, however, remind the reader of what we said on 
p. 340 about “‘things not happening this way.” The economic factor lends 
itself more easily at present to mathematical treatment. That is why we 
treat it here somewhat more in detail, reserving a similar more detailed 
treatment of the other two factors, military and religious, for the future. 
The methodology illustrated by the following is, however, applicable gen- 
erally. 

Let each individual possess on the average an amount s of land. Actual 
holdings will be distributed according to some distribution function. If the 
distribution is due to chance only, as may happen in many cases, then the 
actual holdings will be distributed exponentially (Rashevsky, 1951, chap. 
ix). If F(s)ds denotes the fraction of individuals having a holding between 
s and s + ds, then 


iF (3) ate, (1) 
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Let x denote the minimum amount of food needed per unit time for 
subsistence. Actually x9 varies from individual to individual, but for sim- 
plicity we shall consider x» as constant. Let f represent the average amount 
of food that can be produced per unit time per unit area, when full neces- 
sary work is applied. Then an individual with an amount s of land will 
produce per unit time an amount fs of food. If 


} S<aXee sasor icees (2) 
the individual cannot subsist. On the other hand, if 
Xo 
— 3 
.> 7 (3) 


then the individual may even accumulate a reserve of food. If Xm is the 
maximum amount of food which an individual consumes per unit time 
(x, actually varies from individual to individual), then an accumulation 
will take place, if 


eee 
Seal Seyi 4 
S aegis ae (4) 


Hence a fraction of individuals 


=f F(s)ds= e~*m/sf (5) 


will have a surplus, while a fraction of individuals 
Xo/f 
y= F(s)ds=1— e/a 
s= fF (s)ds=1—6 (6) 
will starve. The remaining fraction 


F,= 1—F,—F3 = e-*m/af + e—*o/ 8s ii) 


of individuals will not starve but will not have an appreciable surplus. 
If No denotes the total number of individuals, then the total surplus 
accumulated per unit time by the fraction F, of individuals is 


E= Nef" (fs—wm)F(s) ds. (8) 


m 


Because of (1), we have 
E= Nof 5 e~*m/fa , (9) 


The total deficiency D of all the NoF; individuals is equal to 


ry /f 
D= Nf (%o— fs) F(s)ds = Noto— Nf3(1— elf) 50. (10) 
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If 
E>D 1) 


then the FN individuals may keep all of the #3» individuals from star- 
vation and still retain a surplus E — D. In exchange for keeping the 
FN 9 individuals from starvation, the FN y individuals may either acquire 
property rights to the land of the FN} individuals (if the notion of the 
land ownership is already established in the society) or put the latter into 
various degrees of bondage or servitude. If (11) is not satisfied, then only a 
part F3N) of the F;N individuals will be kept from starvation. The quan- 
tity F3No is determined by the following considerations. 

It is to the advantage of the FN» individuals to obtain the servitude 
and/or land of as many F3N» individuals as possible for a given expendi- 
ture £ of the surplus. Therefore they will preferentially keep those who are 
less below the starvation point x) than others from starvation. Let xy be 
the food supply of the poorest individuals whom they will help. Then 

E= Ce fslemhds, 


or, because of (9), 
Nof 5e—*m/f2 = Ny{ (%)— x) — f 5) e-*/fi +. f 5 et0/f8} , €12) 
This transcendental equation determines «9. The quantity F is given by 


/ 1 cf —s/s 
R=sf eds. (13) 


The amount S; of land which the 7, individuals will now possess is 
equal to what the F Vo individuals possessed plus that which the FiNo 
individuals possessed originally. Hence 


_ No dae! —s/s No 7 —s/sq 14) 
sao se ds+— ee Sz, ( 


The quantity fs is the average yield & per individual per unit time. If we 
consider that cultivation of extra land was a rather hard task with primi- 
tive means, we may assume that f3 will be very close to %,. If we put 


ff §== 2g, Glsy 
and arbitrarily set x» = 34m, we find 
Pc 37s: F, = 0,03. ; Fz= 0,6); (16) 


that is, about 60% will starve, and contribute to the class of hired individ- 
uals, serfs, or other types of persons with a relatively limited economic free- 
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dom, while 37% of the population will have a surplus. If we put fs = 
X0, Xm = 2X0, we find 


F,=0.13 ; F,=0,5 3+ (Fs = 0.37), (17) 


that is, a much stronger economic middle class. The fraction of individuals 
who receive per unit time more than the amount x of crops is given by 


F(x) = e-*, (18) 


With sufficient land available to the original settling group, & is likely to 
be equal to the average consumption. In that case «/% measures roughly 
how many individuals the individual with “income” x can hire to make 
some work for himself. For «/% = 5, we have F(x) ~ 10-®. For x/# = 10, 
F(x) ~ 10-4. Thus there may be about one individual in a hundred who 
may have a few servants, and only one in ten thousand who may have 
about 10 or more of them. In a small village with a population of a few 
hundred only, there may be a couple of “wealthy” individuals who will 
have 3 to 5 servants, but as a rule none of the wealthier ones. 

‘Instead of hiring men or having them as slaves, for personal service 
only, it is more likely that such a “wealthy” peasant will have the men 
whom he controls manufacture all sorts of luxuries for him, such as pot- 
tery, etc. He may trade some of his food surplus for gems, which may be 
found rather far away from his village. It is known that gem trade has 
been carried over rather large distances in the early peasant village cul- 
tures (Forbes, 1950). This will require the servant or serf of the wealthy 
man to travel to neighboring or even more remote villages, and this travel 
constitutes a medium for the diffusion of information. The fraction of such 
traveling individuals at this stage is also of the order of 10-, since it is not 
likely that the needs of one wealthy individual would require more. Thus 
on the average there will be a couple of traveling persons per village. The 
distance between villages is, as we have seen, of the order of 30-40 km. 
An errand to a neighboring village and back is likely to take a few days— 
say, a week. If the errands are made with equal probability to all nearest 
villages, we may consider the carriers as diffusing particles of information 
(Landahl, 1953) with mean free path A ~ 40 km., and average velocity 
v~ 40 km. week = 2000 km. year. This gives for the diffusion co- 
efficient D ~ 80000 km.? year. The diffusion of a cultural trait—say, the 
use of a particular pottery or ornamentation—will carry that trait in 
t = 100 years to distance 


l= VDtH3X108km. 
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Such estimate may be compared with available archaeological data. The 
estimated speed of diffusion is greater than the one which we estimated for 
a uniformly thinly spread paleolithic population (Rashevsky, 1953), 

The distribution (1) does not arise as such initially. It corresponds, un- 
der certain conditions, to the most probable distribution of a conservative 
quantity among a finite number of individuals. Initially, most likely the 
land may have been divided uniformly, so that each individual had his 
fair share £ of crops, equal to the average consumption. Some individuals 
may accidentally have had hard luck in form of—say, an illness—which 
prevented them from cultivating their lots adequately, or in the form of 
some damage to their crops. They may seek help from a more fortunate 
neighbor, who, as a return for that help, may require either property or 
personal right over the less fortunate. In this fashion a series of random 
accidental factors will gradually lead to either the distribution (1) or to 
some other distribution, for example, the normal one. The general expres- 
sions (5), (6), (7), (8), and (10) hold for any distribution function F(s). 
The method is thus quite generally applicable. 

A different possible mechanism of accumulation of wealth must also be 
considered. This mechanism may have played an important role during 
the early stages of the “metal” culture, when man had just learned to 
produce copper or, later, bronze, and to make tools out of them. The 
village smith, who was the primitive miner and metal worker in one 
(Forbes, 1950), was an expert in a craft not known to all. In those early 
centuries, one smith could easily supply the needs for tools and arms of 
well over a hundred individuals (Forbes, 1950, p. 72). Thus a village of a 
size considered above would on the average have one smith. For illustra- 
tion only, let us make the assumption that each of—say, 200 individuals 
whom the smith supplied—was willing to give him in return 1% of his 
crops. This figure is probably too low in view of existing estimates of the 
price of metal in terms of grain in antiquity. Colin Clark (1951, p. 547) 
finds that as late as the fifth century B.c. 1 kilo of copper was exchanged 
for approximately 40 kilos of wheat. If we allow for one productive agri- 
cultural worker in a family of four, and assume that he made about 1200 
kilos per year, while his need for copper tools was about 3 a kilo per year 
(Forbes, 1950, p. 72), then he would be giving at the above rate of ex- 
change almost 2% of his annual crop. About 5000 B.c. copper was likely 
to be much more expensive in terms of wheat. But even with 1% of the 
farmer’s crop, as assumed above, the smith will receive from 200 persons 
twice as much food as he consumes, and can either barter it for objects of 
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luxury, or use it to hire a personal servant. With a compensation of 3% of 
the farmer’s crop the smith would become a relatively wealthy man in the 
village. 

There is evidence (Forbes, 1950) that in some early communities the 
smith was both socially and economically a prominent individual though, 
curiously enough, in some other societies he was almost a social outcast. 

The formation and persistence of any of the above-mentioned classes 
requires a certain minimal size of a social group. The formation of any 
class, whether purely economic or purely social, requires a definite major- 
ity acquiescence to certain behavior patterns, which frequently may be of 
disadvantage to the majority. This requires in its turn the presence of 
imitative or mass behavior (Rashevsky, 1951a), and that imposes a lower 
limit upon the size of the society, as given by the inequality (1) of chapter 
xili (Rashevsky, 1951) 
a(o+hk) 


Ack (19) 


Ve 

Calculated and observed sizes of communities, which existed at the 
time the first classes were formed, fix approximately the value of the right 
side of (19), and thus may serve as one of the four equations needed to 
determine the four parameters A, a, c, and &. 

3. Formation of cities. We shall define a city as a community in which 
all or by far the largest majority of the population are engaged in non- 
agricultural work. How does a primitive village transform into a city? 

Cities may arise either as political and military, religious, or economic 
centers. Probably all three factors may have been operative simultane- 
ously. For simplicity we shall study each factor separately, and we shall 
begin with the economic factor, for reasons given above. 

We have already touched upon the role of the early village smith. He 
may have been a possible factor in the process of such a transformation. 

Early metallurgy, or any other early craft, could not have arisen in all 
villages simultaneously. It arose most likely in a village which was close 
both to an ore deposit and to a forest, which provided the necessary char- 
coal. Thus while the less favored villages may have been still using flint 
implements, manufactured perhaps by the agricultural worker himself in 
his spare time, the privileged village was supplied with metal tools and 
arms. As we have seen, this was likely to make the smith rich and capable 
of hiring help. He may have decided that with one or two helpers he could 
supply a couple of neighboring villages. If he payed the helpers a fixed 
price, his own surplus would increase by this expansion, and he might ex- 
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pand further, hiring more and more of his fellow villagers, until the whole 
village was hired for a non-agricultural occupation, and supplied a large 
number of other villages. Thus a village is transformed into a small city, 
which numbers only a few hundred inhabitants. The process of expansion 
may continue, and people will be hired from other villages, thus initiating 
a movement from the country to the city. From the earliest stages, when 
only one or two other villages are served, the smith will have to hire men 
not only for mining and manufacturing proper, but also for transport. The 
final size of the city will be reached when transport difficulties become so 
great that a further expansion is economically unprofitable. 

While, as we said above, in some primitive communities the smith 
might have been a socially and economically privileged individual, in 
others he might have been almost a social outcast (Forbes, 1950, p. 63) 
and would have no chance to start or operate the above-outlined process of 
expansion. In this case the process may have been carried out by the 
“wealthy” man of the village, who could hire the smith. The mechanism 
remains essentially the same. 

We shall now proceed to a mathematical description of the above- 
outlined process of formation and growth of a city. We shall speak neither 
of the smith nor of the ‘“‘wealthy”’ man, but of the “enterpriser,’’ who may 
be either of them. 

Let 1/p’ denote the number of individuals that can be supplied by one 
primitive metallurgical worker. As we have seen, 1/p’ is of the order of 
100, so that 

p'~107-. (20) 
As long as the village is located at the source of the ore and near to a sup- 
ply of charcoal, the problem of transport of raw material does not arise. 
As the enterprise expands, more raw material is needed, and this may have 
to be carried over a considerable distance by burden-carrying animals or 
by boat. We have mentioned so far two kinds of raw material, ore and 
charcoal (or wood). Clay is also necessary to build the furnace. In general 
let us consider any number of raw materials. Let c; be the amount of the 
ith material required per unit time for one person actually engaged in 
manufacturing and /; be the distance over which it must be transported. Let 
there be V), persons engaged in active manufacturing and let N; be thenum- 
ber of individuals engaged in the transport of raw materials. If 2; is the 
average speed of transport, which may vary for different raw materials, 
then the number of round trips per unit time for the 7th item 1s 0;/2l;. 
If v is the load per attendant (Rashevsky, 1953), then the total number 
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N;, of attendants needed for the ith item is 


yi a 2 Neal (21) 
: Vid; 
We have Le 
Me= N= Mey (22) 


Let W be the total number of individuals supplied, or the total number of 
consumers of the finished product. The total number of individuals needed 
in the process of manufacturing, including transport of raw materials, 
is V,, + N;. Hence the ratio 1/p, = N/(N, + N/) denotes the number of 
individuals supplied per individual engaged in production. From (22) we 
have, remembering that p, = V,/N 


Ne NO. 2¢,l; 
meet AN tic BY) 23 
ees ee CED “= (23) 


Thus with the need of transport of raw materials, the number of individ- 
uals in the producing industry per consumer of finished products increases. 

Put 
N,=N'+Ni. (24) 


This is the total number needed for “production.” 

In order to supply WV customers, the enterpriser will need VV, hired men 
for production, and a number JN; of hired men for distribution. Let the 
transport speed in distribution be v and the load per attendant v, and let 
the same values hold for the transport of rural agricultural products to the 
city. Put 

vu=4u, (25) 
and let us first consider the land transport only. Let the V consumers in- 
habit an area of radius r. Consider also the oversimplified situation in 
which the same transport individual delivers the finished manufactured 
product to the rural population and brings back to the city the agricul- 
tural products. Roughly such an individual will make v/2r round trips per 
unit time, and therefore the total amount of agricultural products brought 
to the city per unit time is 


N,u 
Pa (26) 


If every rural individual is supplied with the manufactured product, then, 
if D is the density of population, we have 


N=rr'D, (27) 
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and since 
Ng = pal ; (28) 
therefore, equating (27) to (28) and solving for r, we find 
Ny 


r : 
T pu D 


(29) 


Let each rural individual produce on the average an amount #, of food 
products per unit time and let him be willing to give the fraction ¢ of it for 
the needed supply of metal products. The total amount of food trans- 
ported to the city is then 


Ny€?Pr 


; (30) 
Pu 


and this is equal to (26). Hence, equating (26) and (30) and introducing 
(29), we find 


(31) 


Of the total amount NV ,e€,/p. of agricultural goods received per unit time, 
the enterpriser must use the amount J ,c, for feeding the V, manufactur- 
ing workers, and NV ,c, to feed the distribution workers, if c, is the per 
capita food consumption per unit time. Hence his surplus G is equal to 


GaNvhe_ (N,-+N,) Ce. (32) 


uu 


Introducing (31) into (32), and putting 


fr _ ¢ = a(>0) (33) 
Pu 
2€p; Cr 
Web aan Gs 
we find 
GaN NY (35) 
U Pp 


The enterpriser will hire such a number of men N > as to maximize G. 
But G has a maximum for 


= (2a8y. (36) 


If he hires N} for production, he will hire a patties Nj for distribution. 
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That number is given by (31), into which we introduce (36) for V,. Hence 


8 aru? 
ee S71 
NS=97 Bic,” (37) 
The function p of the total population, which is engaged in distribution, 
is 
* 

p= N; =e 2 a pu ; (3 8) 

No 2ap,+3c, (1+ p,) 

mompete we) tie es 
Pu 


If we include all transport workers as belonging to the city population, 
then that population is equal, because of (31) and (37), to 
12a?2c,+8a3 


2 3 
27 Bc, u? . (39) 


n=N*+N¥= 
As in the simpler case treated previously (Rashevsky, 1953), 2 increases as 
u? = (v)?, 

In this oversimplified picture we considered only one kind of “indus- 
try,” the making of metal. If several different types of industries develop. 
the result, 2n the first approximation, will be the addition of numbers of all 
individuals involved, and a corresponding increase in the size n of the city. 

If a city is located on a river or on a sea shore, then this has the effect of 
substantially increasing wu and, therefore, such cities will be larger than 
inland cities. If the city supplies a population, part of which can be 
reached by land, part by water, then the value of u is some appropriate 
average of the values u and mu, for land and sea transport. This average 
depends on the geometry of the region supplied by the city. 

Of special interest for future studies is the fraction p Of the total popula- 
tion which is engaged in distribution. This fraction may be taken as a 
measure of the dissemination of the information of any kind, in particular 
of information about nonconformists and their new ideas. For the fraction 
p essentially represents the fraction of the total population that travel 
from place to place. In ancient times, in absence of such means of mass 
communication as press or radio and with the use of a postal system con- 
fined to a very few individuals, any information about what was going on 
in other cities or villages could be obtained almost only through traveling 
individuals, who were engaged in distribution. If p = 0, then every indi- 
vidual is confined to his immediate circle of acquaintances, and the flow of 
information in the country is at a minimum. Information may still pro- 
ceed through occasional visits amongst neighboring communities. Its 
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spread is then described by a diffusion equation, similar to one discussed 
above. [Also Rashevsky, 1953, eq. (88).] This process is very slow. If 
p = 1, then we have the maximum possible spread of information, in the 
absence of methods of mass communication. As we have seen (Rashevsky, 
1953), the rate of spread of information is an important factor in deter- 
mining the rate of historical development. Hence the importance of the 
quantity p. 

A glance at expression (38) show that p does not depend on w, in other 
words on v, though both NV; and N/ vary like w?. It may be thought that 
this result is tied to the special assumption expressed in equation (27). 
This expression assumes that the city is located in the center of a circular 
region. If we consider a linear arrangement, like a very narrow strip of 
land along a shore or in a narrow valley, then N varies not like 7?, but like 
r. In general we may put NV « r2,where1 < q < 2. If, however, we repeat 
the whole argument with this assumption, we still find that p is independ- 
ent of wu. This result is rather general. 

Consider, however, a city located on the shore of the sea and selling in- 
dustrial goods to foreign countries across the sea, in addition to serving a 
region of radius r in its own country. To serve that region, the city will 
need a number NV} of manufacturing workers, given by (36), and a number 
Nj of transport workers, given by (37), and will contribute the value of p 
as given by (38). To serve the countries overseas, it will need a number 
N3’ of manufacturing workers, where Nj’ is given by the expression (36), 
in which we substitute for uw the value u, for sea transport, and for D (in 
8) the value D’ as referring to average density of foreign population for 
the total area covered by sea travel. Thus 


2au 
a 40 
Ny ) (40) 
where 
2p, Cy (41) 
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More generally we could also assume different values for e, ?,, and Gr over- 
seas. We shall, however, for the present omit this further generalization. 

Similarly, the number N7" of individuals engaged in sea transport 1s 
now given by [see eq. (37)] 


8a3uz (42) 
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However, while in the case of inland cities the number JV of transport 
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individuals was recruited from those who inhabit the area served by the 
city, in the case of overseas trade, the number N*’ is recruited from the 
individuals inhabiting not the overseas area served, but the home area. 
Hence now the total fraction p’ of individuals engaged in transport is 
given not by (38) but by 


yr* yr’ 
p’ —_ ae + N; (43) 
+NF +N? 
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where N;, N35, and N;" are given by (37), (36), and (42) respectively. 
Hence 
a oe? (44) 
If we consider a city on the borderline of the country, where that border- 
line is on land, we can make a similar argument, and introduce the number 
N7"’ of individuals engaged in distribution to the outside or foreign coun- 
try, but recruited from the population of the country considered. In this 
case, however, N;"’ is given by the same expression as N7, that is, by 
(37), because land transport and not sea transport is involved. Because of 
a possible difference in the population density D of the foreign country, 
the value 6” of 8 may be different in V7” than in N*. The corresponding 
value p” will be greater or less than p’, depending on whether N*”’ is 
greater or less than N?’. But, other parameters assumed to be the same, 
we have 


2 2 
, ae ©. Us 
Nay alec ait ee (45) 
Since 8 varies as D-/?, therefore inequality (45) is equivalent to 
D'u2> D''u2. (46) 


As we have seen (Rashevsky, 1953), u2 ~ 1042, Hence, unless D’” ~ 
10* D’, inequality (46) is satisfied, and p’ > p’’. 

If we have a shoreline along an ocean, where the nearest countries 
across the ocean may be out of reach of the available navigation methods, 
then D’ is almost zero, and p’ = p. The effect of the shoreline is negligible. 
If, however, as in the Mediterranean, a large number of countries is ac- 
cessible for trade, p’ will be appreciable. A glance at the map shows that if 
we take the population around the Mediterranean, in a strip approxi- 
mately 100 km. deep, and estimate the population density by including 
also the area of the sea, this gives us a value of D’ between 0.1 D’” to 
0.01 D’’. Thus p’ can be many times as large as p”. 
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Now consider a country of total area S, and a sea shoreline £5. The aver- 
age distance between cities is of the order of 2r*, where r* is obtained by 
putting V, = N; in (29). Hence, because of (36) 


ee fe 
TEV ra: Ga 


The number of seashore cities is then approximately 


ES 
mM; = ors . (48) 
The total area S, served then is approximately £57*, since it is represented 
by a strip of length &S, and depth 7*. Therefore the area left is 


S,=S—S,=S(1—ér*). (49) 


The fraction of transport people in the area S; is equal to p, that in the 
area S, is equal to p’ > p. Hence the average value of p for the whole 
country is 


* of aS: 
Se ahd 8 Fh ga (50) 


and increases with £. 

If, like in Greece or Italy, the distance from the deepest inland point to 
the sea is less than 7", the city will arise on the seashore, since this maxi- 
mizes G. There will be no significant inland cities. For such countries p 
will be much larger than for countries like India or China. 

We shall now briefly discuss the problem of growth of a city from a 
small village to its maximum size given by (39). As long as V, + N;zis less 
than the population n, of an agricultural village, there are no cities prop- 
erly speaking. A city “is born” at the moment when V,+ N, = n,(~ 
500), that is when all the population of the village is hired for production 
and distribution. From that moment, by our definition, the village is trans- 
formed into the smallest city. 

An increase in V, and JV, can, as we have seen, occur by the enterpriser 
hiring new men, using a part orall of his gain G. In the model adopted here, 
the gain G is measured in terms of agricultural food products, and is given 
by (35). Since these products are collected normally once a year, therefore 
the maximum yearly increase of V, + N, is obtained by investing all of G 
in hiring new men. Since the consumption per individual is c,, therefore 
we find, expressing / in years, 


d(N,+ Mi) _G (51) 
dt Ce 
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Introducing (31) and (35) into (51), and remembering (34), we find 


d (x, +2") 4 B 
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Integrated with the initial condition NV, = 1 for ¢ = 0, equation (53) 
gives 
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Once we know J, as a function of ¢, we know also NV; as a function of #, be 
cause of (31). The time at which the village is transformed into a city, 
counting from the first smith establishing himself professionally in a vil- 
lage, is given by the root of the equation 


N,p(t) + Ni (t) =n, (55) 


solved for ¢. The time at which the city reaches its maximum size is ob- 
tained by solving for ¢ the equation 


N,(t) = Np , (56) 
where NV; is given by (36). 

In principle it should be possible to obtain from historical data, in par- 
ticular from some records, information about the values of p (fraction of 
population engaged in distribution), or about the time it took certain cities 
to reach a certain size. [Cf. rough estimates in Childe (1952), p. 94.] As 
mentioned before (Rashevsky, 1953), we thus may be able to ““postdict”’ 
some historical data. 

A discussion of plausible numerical values of the different parameters in 
equations (38), (39), (50), and (54), and of the numerical results to which 
those equations lead will be postponed to a subsequent paper, after we 
derive similar equations for cities which are formed as military-adminis- 
trative or religious centers. 

In the above discussion we have made several oversimplifications, 
which in the further development of the theory must be gradually dis- 
carded. One important oversimplification is the assumption of a single dis- 
tributing center—the city itself. Conditions of transport may be such that 
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it is easier and cheaper to establish a number of smaller distribution cen- 
ters, which form smaller cities. Those may exchange goods directly with 
the villages. This introduces a specific structure into the transportation 
net, which we hitherto considered as amorphous. A connection with the 
recent work of Martin Beckmann (1952) is thus established. It is also 
likely to lead to expressions for the distribution of city sizes. 

This problem, as well as the theory of cities which are formed as ad- 
ministrative or religious centers, will be the subject of future communica- 
tion. 


The author is indebted to Drs. George Karreman and Anatol Rapoport 
for a careful reading and critical discussion of this paper. 
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AN AGE-DEPENDENT STOCHASTIC MODEL 
OF POPULATION GROWTH 
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A stochastic model of population growth is treated using the Bellman-Harris theory of age- 
dependent stochastic branching processes. The probability distribution for the population size 
at any time and the expectation are obtained when it is assumed that there is probability 
(1 — ¢),0 S o <1, of the organism dividing into two at the end of its lifetime, and probability 
o that division will not take place. 


The theory of stochastic branching processes is concerned with the 
study of processes representing the evolution or development of systems 
whose components can reproduce, be transformed, and die; the transitions 
being governed by probability laws. W. Feller (1939) was the first to 
apply these methods to biological problems in a paper on stochastic models 
of population growth. Since that time the theory of branching processes 
has found many useful applications in biology (cf. Kendall, 1949). 

Recently R. Bellman and T. E. Harris (1952) have developed a theory 
of age-dependent branching processes which appears to have important 
applications in biology and in physics. In this theory the generation time, 
7, from the birth of an organism until its own subdivision is a random 
variable with the general distribution G(r), 0 <7 < ~. At the end of its 
own life the organism is transformed into ” organisms with probability 
Qn) % > 0, each having the same distribution of generation times G(r). 
The Bellman-Harris process is formulated as follows: let X (¢) bea random 
variable representing the number of organisms at time #, and define 

p,(t) = Prob[X (¢#) = 1], r20 
and 


$(s, 1) = >> p-(0) s%,| 5) <1 (1) 
r=0 


as the generating function for the number of organisms starting with one 
at t = 0; [(s, #)|" is the generating function for the probabilities of the 
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number of organisms at time ¢ starting with n at ¢ = 0 and assuming inde- 
pendence of the organisms. Bellman and Harris have shown that @(s, #) 
satisfies the functional equation 


o(s,1) = [blo (st—7)1dG(r)+51-G)I, (2) 
0 
where 
h(s) = >) as". (3) 
n=0 


Assuming G(t) to have a density function of bounded total variation, 
t foo} 
G ti) = d d <o, 
=f g(nar, fi \date)| 
we can write (2) as 
t 
o(s,) = f hio(s,t—7)) g(r) dr +s11-GO]. 
0 


Differentiation of (2) with respect to s yields integral equations of the 
renewal type for the moments, the properties of which can be studied 
using well-known methods (Feller, 1941). 

Bellman and Harris in their studies have restricted themselves to the 
study of a pure birth process (i.e., a process in which the random variable 
representing the population size can only increase whenever it changes) 
in which transformation or division is always binary. In this note we con- 
sider the case in which at the end of the life of an individual there is prob- 
ability o, (0 < o < 1) that division will not take place, and probability 
(1 — a) that division will be binary. Hence we have for (3) 


h(s) =o+ (1—¢) s?. (4) 


If we assume G(t) = 1 — e~*, where X is the expected birth rate per or- 
ganism per unit of time, the process is of the Markov type and (2) be- 
comes 


News -f' (1-woj dis, b= chol he ord + eek eee 


With G(t) as defined above the state of the system at any time ¢ depends 
only upon the organisms existing at that time and is independent of their 
ages. The probability that an organism existing at time ¢ will be trans- 
formed between ¢ and ¢ + At is \At + o(AZ), independently of age and 
absolute time. We could assume some other form for G(#), which would, 
perhaps, be more realistic in biological applications; however, the ad- 


— 
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vantage of the above form is that it enables us to obtain a solution to (2) 
with little difficulty. The case where G(t) is a convolution of & distributions 
of the form 1 — e~t has been treated by D. G. Kendall (1948). Here the 
organism passes through & stages before it is transformed. 

To solve the non-linear integral equation (5) we write £ = 1 — 7, so 
that 


t 
eg ( 5, 8) zip {(1—¢) $2(s, —£) +o}vAede+ s. 
0 
Differentiating the above with respect to ¢ we obtain the differential equa- 
tion 


d¢ 
Ot 


The solution of (6) is 


=A{ (1-9) ¢?—¢+>¢}. (6) 


2(1—0c)¢-—1-— 


Typ 8 HIST KS), (7) 


where we have put p = (o? — o + 4)”, and K(s) is to be determined 

using the initial condition ¢(s,0) = s. After determining K(s), and solving 
(7) for ¢, we have 

2(1— 0) s—14+20], n,, 

Pathe Mato) Se a) rao ae : 


SRT TS eres ae 


(8) 
2(1—o) }1 Fis) s—1.—= 2p 


To obtain the probability distribution for the population size at any 
time ¢ we must expand (8) in powers of s. Let us rewrite (8) as 


at+a’s 
Sigs Tale eae are (9) 

where 

o = 4 p?— 1+ (1 — 4p?) e—"), 

a’=2(1—c) [1—2p— (1+2p) e~”*‘]. 

p=2(1—s) [(1—2p) Wipe Ata Meas) 

=4(1—¢) (o—1) [1—- Ee rpt] , 

Now 


¢ (5, t) =f ata’s)-(1-$s) =5) (ata’s)- >) st 


therefore 


p= ACG), Bt. (10) 
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In (10) we have put A = (a/B + ap ye 
We now consider the case when ¢ = 3. When a has the above value the 
solution of (6) does not hold since p = 0. In this case we have for equa- 


tion (6) 


Op _ Xr 
ayer Be Sid (11) 
hence , ; iy 
1 vAE— 1+ (2—At) s 
ae )-Gant aadi( ae Ys (12) 
1+ Xt 


Proceeding as before we have 


POS ret +==*) (Gx): (183 


As stated earlier, the expectation of X(¢), or the mean population size, 
is obtained by differentiating (2) with respect to s, and evaluating the 
derivative when s = 1. Putting EX(#) = y(t), we have 


vo =(F)_ f¥e- dG) + 1-G01. 


Therefore 
t 
v(t) =2(1—0) f ¥(t—2) ded te, (14) 
0 
This integral equation reduces to the differential equation 
d 
Ho y(1-20)¥, (15) 
whose solution is 
v(t) = EX (t) =X (0). AC. (16) 


In (16) X(0) is the initial population size. It is clear from (16) that the 
mean population size will depend on the value of o. For « > 4 it decreases 
as ¢ increases, for a < 4 it increases as ¢ increases, and for o = 3 it has the 


value X(0), hence the mean population size is stationary. As ¢ approaches 
infinity we have 


and 


lim EX(!)=0, o>4 
(17) 


limEX(t)@, o<}. 


We remark in closing that upon putting ¢ = 0 the process treated 
above becomes the well-known Furry process, equivalent to the case 
treated by Bellman and Harris, which is of interest in the theory of cosmic 
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ray showers (cf. N. Arley, 1948, p. 93). The Furry process has recently 


been used by S. Iversen and N. Arley (1950) as a model of the growth 
process of tumors. 
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The transmission of some information or behavior pattern is treated as a flow of “‘particles’’ 
which execute random motions over a population of individuals and which may multiply or dis- 
appear. Equations are derived for the number density of these ‘‘particles” and from this is cal- 
culated the number of individuals through which the ‘‘particles” have passed. The results are 
applied to a number of situations such as 1) uniform spatial distribution with multiplication 
factor decreasing with time because of loss of interest or confusion of the information, 2) mul- 
tiplication factor constant but the rate of spread decreasing with multiple hearings, 3) one- 
dimensional region with a small starting region with or without an absorbing barrier, 4) two- 
dimensional region with absorbing barrier, 5) continuous sources of information within a small 
region in one dimension, 6) uniform spatial distribution in which individuals do not respond to 
more than one hearing. 


We shall consider the following model for the spread of a state (idea, 
piece of information, a behavior pattern, communicable disease, the ac- 
tivity of neurons, etc.) over a population. Let “particles” diffuse through 
a space occupied by a population. The diffusion will follow the well-known 
laws of the diffusion of gases, except that the particles can also multiply in 
addition to diffusing. In fact, a particle multiplies each time it hits an indi- 
vidual of the population, and the multiplication factor F may be greater 
than, equal to, or less than unity. 

Translated in terms of information spread, this process implies that an 
individual transmits a piece of information each time he hears it, and F is 
the average number of times he transmits it. In general F can be a function 
of time. In particular, if the individual transmits the information upon 
first hearing only, our F is equivalent to the parameter a in the papers on 
neural nets and contagion spread of R. Solomonoff and A. Rapoport 
(1951), G. H. Landau (1952), and Landau and Rapoport (1953). 

There are, of course, other interpretations of our model. For example, 
the particles can be regarded as parasites which multiply within blood 
cells and diffuse when the last cell bursts. On the other hand, they may be 
regarded as impulses traveling in a random neural net. 
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The method pursued here is a generalization of the method used by 
Landau and Rapoport in the sense that repeated transmissions by the 
same individual are allowed (the simple transmission, i.e., the case with 
‘perfect immunity” being a special case) and spatial factors are also 
treated. However, our results are, in a way, more restricted than those of 
Landau and Rapoport in that the “private time,” i.e., the time between 
hearings and tellings of each individual is not considered. 

It may seem at first that our model violated an important principle in 
the spread of a state. Usually in a spread of a state, the “state” is not 
thought of as necessarily “leaving”’ one individual when it is acquired by 
another, whereas our particles, behaving like material particles, do so. 
It must be emphasized, however, that the particles serve only as ‘“‘coun- 
ters’ in our process. It is not where they are but where they have been 
which indicates the number of individuals which have been involved. On 
the other hand, the number of “particles” is a measure of the intensity of 
activity in the population at any given time. 

This model of the spread of a state makes it possible to avoid the major 
difficulty which immediately arises in the usual treatments, i.e., the non- 
linearity of the equation for the rate of spread, which is due to the fact 
that only that part of the population which the state has not yet reached is 
available for recruitment into the affected population. However, we can- 
not simply count the number of individuals who have been reached by the 
particles [n(P, ¢) in (3) below], because this would mean that any individ- 
ual to whom—-say, in the case of information spread—the information had 
been told more than once would be counted more than once. A correction 
for this error is made by using (4) below to give the fraction of the popula- 
tion who have actually heard the information at least once. 

Consider the transmission of an idea, some information or a behavior 
pattern from one individual, 7, to another, 7 + 1, as a “particle’’ which 
moves randomly in a space occupied by Yt individuals. Let 7 be the aver- 
age interval between the time that the information is heard by a given 
individual 7 and the time that it is repeated by i to i + 1. Let X be the 
mean value of the distances /, where / is the distance between the position 
in space at which the information was heard and the position at which it 
was heard by the next individual. We shall suppose that the individuals 
are relatively stationary except perhaps when transmitting information. 

Let F be the mean value of the ratio of the number of repetitions to 
others to each hearing, e.g., if each tells two others on the average F = 2. 
Let m(P, t) be the number of particles occurring per unit of space (area or 
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length) about a point P at time ¢. If F = 1, the space integral of m is inde- 
pendent of time. If the particles move randomly over a one- or two- 
dimensional (vy = 1, 2) space, we may introduce a coefficient D which will 
determine the rate of spread of the particles in time and space 


LiSieetacens (1) 


The last expression in (1) is exact if the distances / are distributed ex- 
ponentially. Equating the rate of increase of m with time to the flow into a 
given region due to diffusion plus the increase when F > 1 or decrease of 
F < 1, we find 
om = DVM to) (2) 

In writing (2) we have assumed that D is constant, that the average 
value of the vector distances / is zero, and that even if particles meet on 
the same individual, each particle multiplies unchanged. The first assump- 
tion is satisfied if the mean distances that the particles jump, as well as the 
mean of the times between jumps, are independent of position and time. 
If the assumption does not hold, the quantity D should occur under the 
operator (Patlak, 1953). If the second assumption does not hold, one must 
subtract the divergence of the product of m/v7 and the mean value of the 
vector distances /. The effect of the third assumption is not important 
unless m is fairly large. This point will be discussed further in Case 2. 

For the moment we shall suppose that (2), together with sufficient 
boundary conditions, is solved for m(P, t). Consider a given region about 
the point P. From the time ¢ = 0, prior to which m(P, ¢) = 0, up to the 
time ¢, the total number (P, #) of hearings per unit space is given by the 
sum of m(P, t) over units of 7, or 


n (P, t) =~ fom (P, 1’) dl’. (3) 


In this expression, as also in (2), we assume that the particles are not 
synchronized in their movements; otherwise, for example, one would use a 
sum instead of an integral in (3). If one synchronized the starting of “par- 
ticles,” there is then a short period of adjustment to randomness. The 
effect is chiefly that of a delay by about 7 together with a spreading out in 
time of the discontinuity. This can be illustrated for simple cases. Let 
D=0 and F = 0, ie., no repetitions. Then the solution of (2) can be 
written, since m and are independent of P in this case, 


m(t) = moe" n(t) = mo (1— e-"*), 
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so that it is as if the mp particles were introduced between the individuals 
but do not contact them until 7 later, on the average, the contacts being 
spread out over a few units of 7. As another illustration let D = 0 and 
let F = 1 from ¢ = 0 to ¢ = 7, then F = O thereafter. In this case if 
synchronization was maintained, then m(0) = mo, m(r) = mo, m(t > 7) =0, 
n(0) = mo, n(t > 0) = 2m. However, from (2) and (3) we find m = mo 
and # = mot/r fort < 7, 


m=mee—/", n= m,(2— ee—*/7) for t=7. 


Here again the discontinuities in m(#) are smoothed out and there is a 
delay of about r. 

If in the same region about P there are NV (P) individuals, then the num- 
ber V(P, t) of these V7(P) who have been involved in the transfer of in- 
formation will, for small 2(P, ¢), be nearly equal to 1(P, t) but for larger 
values it will be less than u(P, #) because of the duplication. Let the prob- 
ability that any one individual has been involved in r events up to time ¢ 
be P,(/N 7). If this can be approximated by a Poisson distribution 


e-wNT (s:) 
Pp (+) = AN tf 
- Ne r! ; 


then the fraction N/N, of individuals who have been involved in the 
transfer process, which is one minus the fraction not involved, is given by 


eS ig la pealy ad ladle a pl Rm 
an Psa) =! og: (4) 


If Nr is not large enough for the Poisson distribution to hold, then we have 
the following instead of (4) 


N 
Speen i es 

For simplicity we have assumed a uniform spatial distribution of individ- 
uals. If the individuals, after a time #, tend to forget the fact that they 
were involved in the transfer of the information in question and if the 
function f(t), monotonically decreasing with #, is a measure of the amount 
retained after a time /, then (3) can be modified by introducing f(é) under 
the integral to obtain the number which still remember the information. 

Case 1. Uniform spatial distribution: F = Fye-*t. Let the events be 
randomly scattered over space and let m(#) be the space integral of m(P, t). 
In the spread of some types of information, the frequency of repetition, 
1/7, may decrease with the time of initiation of the process. Also the mean 
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value of F may decrease if the information is dated and the interest in the 
information decreases. Since the results are much more sensitive to the 


change in F than in r, we will consider only the former. Thus if we have, 
for example, 
F= Fy c= ? (S) 


and since m is independent of P in this case, expression (2) becomes 


dm(t) _ (Fye~**— 1) m(t) 
ae Sage ae me Ko) 


Thus, if mo = m(0), 8 = 1/ar, then from (6) we find 
m (t) = my el FBU-e-%)—t/7] | CT) 


Introducing expression (7) into (3), we obtain [n(¢) = n(P, £)] 


m,e* 8 
at 


n(t) = Fup) ef ey Ady. (8) 


If we simplify matters somewhat by letting F = F)(1 — a) and stop 
the integration at ¢ = 1/a when F = 0, then instead of (8), we find for 
t < 1/a, if y? = 1/aPo, 


Meo Z (1—F,)y+t/yr : 
n(t) = e(Fo—1) prarye fi ev /2d 9 
V aF yt ‘Giese F,)Y ", y ( ) 


N(t) then being given by substituting 2(¢) into (4). After a time ¢ = 1/a,” 
will increase by mo, most of the increase being in the next interval r. In 
either case for a short time (¢) is small and hence W(t) is approximately 
equal to 2(t). Therefore d’N/dt? at t = 0 equals d’n/dt? = dm/dt. If 
dm/dt > 0 at t = 0, or if Fo > 1, W(é) is concave upward, otherwise not. 
If a= 0, N(~) = Nr unless Py < 1. 1f Fo < 1 
N (@) —m,/N 7A F,) 
Ve =1l-—e z : (10) 
This result may be also used for the situation in which the process does 
not change but the information is distorted with successive repetitions in 
such a manner that for each t/7 the probability of being lost is p;. If pi is 
small the probability of information not being lost at ¢ is given by 


e7Pit/t P 


Hence, if a = f,/7, the desired solution for this situation is given by ex- 


pression (8) or (9). ad 
Case 2. Information spread with uniform spatial distribution. Rate of 
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spread decreases with multiple hearings. Constant F. If an individual, after 
repeating the information F times, on the average, does not respond to 
subsequent hearings for a time 7, then the factor F is effectively reduced 
by a factor 1 — N’/Nz, where N’ is the number which have heard the 
information within time 7. But ”’, corresponding to N’, is the integral of 
m/r from t — T to t. For not too large T this is approximately Tm/r. 
Thus if we introduce a = (F — 1)/r and 8B = TF/(F — 1), andif Fisa 
constant, we find instead of (2) 


am) = am(1—Bm), (11) 

so that 
™ ) = EOE By = Baie @2) 

Hence 
n (t) = tag bs (Bm + e-*'— Bm e*) (13) 


and V(t) is given from (13) and (4). In this case d?V/dt? > 0 only if 
a > OorF > 1. Butunlessa < 0,N(©) = Nr. If F < 1, ora < 0, then, 
if Bm) <1, anda = —a 


N (@) ={— (f=<ipien) er ee (14) 
Nr 


Thus the limiting fraction exposed to the information increases with F 
and with m/Nz ~ No/Nr. Unless Bmp is an appreciable fraction of one, 
the effect of 8 is negligible. 

Case 3. Finite or infinite one dimensional region. Approximation solu- 
tion. Initial distribution: mo ts finite only near the center of region with cen- 
tral symmetry or small region near a reflecting barrier. Let the initial number 
of events be confined to a region near the center of a region p = 0, where p 
is distance from origin, in such a manner that m(p, t) = 0 for p > rp and 
m/(p, ¢) is linear in p for p < ry) and decreasing with p, being zero at p = fo. 
Let m(0, 0) be denoted by mo. If M(é) is the integral of m(p, #) over all p, 
then M(0) = M® = moro for a one-dimensional case 7 = 1 or M® 
= 7mor,/3 for a two-dimensional case j = 2. By the approximation method 
(Landahl, 1952) the solution of (2) is given by [a = (F — 1)/7] 


at p 
— ite Yee 
a ee ae ee ee oe ee ee (15) 


8 DiNi/2 
ister) 
0 


r 
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for p < Vr? + 8Di, but m(p, t) = 0 for larger values of p. Thus we find, 
if p = p/ro, t = at, and A = 72/8D 


E i my Ai/2e—A At+t 
CE eed (EUR pArng-itin) de. (16) 


QT Ap? 


The integrals are tabled functions, e.g., the exponential integral, error 
function, etc. For example in a one-dimensional case and ifa = —a < 0, 
A = —A, we have 


Pav _ paAtat a 
n (p, t) Shes [ v2 fj —ewdy — soe ny — dn| (17) 
p 


Va/4D Va/4D 


or, if Z(*) is the integral of the normal distribution from — © to x and 
—Ei(—x) is defined for positive « by the integral from x to ~ of e-"/n 
with respect to 7, then for p < +/r2 + 8Dt 


n(p, t) Sule [avez {za +au) _2(exa)y 


=) ee 


n(p,t) =0 — for p=Vri+8Di. 


_eva \Ei(— A —at) —Ei(- 


After an infinite time we have 


n(p, 2) = mM V/A eA [2vz} jrez(ome)) 


x ay (19) 
pVvA Va 
nF ro D IK 
so that 
n (0, @) = VrA ed (20) 


and m approaches zero fairly rapidly for p appreciably greater than 
./4D/a. Hence the “radius of action’’ increases with the diffusion co- 
efficient but decreases with a. 

When a > 0, V(p, #) ultimately becomes one everywhere, the function 
having non-zero values at any p whenever /r} + D8 is greater than the 
chosen p. ; 

If the term Bm? is added to this differential equation the effect is small 
ifa < 1. Butifa > 1, the effect is to produce a constant velocity for the 
wave front as it progresses, after certain initial processes are over. 
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If the source is at a reflecting barrier, then replace My by Mo/2 and let 
the function start at mo, dropping linearly to zero at ro. The equations are 
then the same as above. The reflecting barrier corresponds to the situation 
in which an individual at the barrier transmits information to the same 
degree as any other individual within the region but always back into the 
region. 

If the region is finite and of width 2A with reflecting barriers at + A and 
initially m(p, t) is symmetric about the center of the region, the above 
equations for m hold approximately until ¢ = ¢*, where r§ + 8Di*? = A?. 
Thereafter the linear approximation becomes even simpler (cf. Landahl, 
1952, Case 2). If the barrier absorbs information (cf. next section), then at 
Am(p, t) approaches zero and the problem can be solved by integrating 
the solution for the analogous case (Landahl, 1952, Case 2). 

Case 4. Absorbing barrier at +A, F = 1. One dimension. Uniform initial 
distribution. Consider next a one-dimensional region of length 2A in 
which, at time ¢ = 0, mp» “particles” are initiated uniformly over the re- 
gion. Let the individuals near the ends of the region make fewer contacts 
because of the absence of individuals on one side. We may suppose that 
there are fictitious individuals on the outer side who absorb the “par- 
ticles,” i.e., they listen to the information but do not repeat it. Thus, on 
either side of the region the value of m will be zero at some effective dis- 
tance of the order of magnitude of \. We shall ignore this effect because 
the region is large compared to \ or suppose some suitable correction has 
been applied in order to use the following as a first approximation. 

Let m/(p, t) = mo(1 — p/r) where p = |x|, « = 0 at the center and 
where 7 is a distance which depends upon time. Let F = 1 in the differen- 
tial equation (2). Then at any time ¢ < ¢* = A?/4D we may approximate 
m/(p, t) by 


m(p,t) = my, OSpSA-r, 
= 21 
m (p, t) po Bs A= ssn (21) 


Since the number of “particles” which have left the region equals the 
integral of the rate of leaving, we have 


_ f* Dmdt' 
$7 (t) my = i waa (22) 


so that on solving for r we find 


r=2/Di. (23) 
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Hence, for A — 2\/Di < p < A, 


ns A—p A? 
m (p, t) OORT DI ‘ap? (24) 
For ¢ > ¢* let m(p, t) be given by 
= p 
m (p, ) =m*(1-2), (25) 
Then, instead of (22), we find 
t * 
mA —}m*a= f, am dt (26) 
Solving for m* we find 
m* = My e~2Dt/ A*+1/2 (27) 
Hence from (25) for ¢ > ¢* 
m (p, t) = my el en2D1/ a (1 -£). (28) 


Therefore integrating m/(p, t’) over t’ from zero to ¢ we find n(p, t) to be 
given by 


t tS (A— p)? 

m(p,t) = mo, EES 
my(A—p) Vi mo(A— p)? (Am pe A? 
,t) =e Et i eh RI SY i < ool 
AGM 7VD ingest Deg =. 2D; 


= Mo “ie _ My (A?— Ap) / 2Dt 
ee Senate) (AS Mica Pe 


If (0, ~) is sufficiently less than Nz and if mo = mo/Nr, Ni = 2ANz, 
1 = 1/pN7, p being the probability of an effective contact, then from 
(1), (4), and (29) we find the fraction of all individuals within the region 
who have been exposed to the information is given by 


5 
0 pm Ni. 


If the region is long strip of width W and length 2A, the long edges being 
reflecting barriers, the corresponding expression is the same as the above 
if p is replaced by pw/W, w being the effective width of an individual with 


respect to encounter with other individuals. 
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Case 5. Absorbing barrier at p = A. Two dimensions. Uniform initial 
distribution: F = 1. We consider next the same case as that above except 
that the region is two-dimensional. Let m(p, 0) = mo(1 — p/r) and let 
F = 1. Then let m/(p, ¢) be defined again by (21). Instead of (22) we have 


A — 
Btad Ace ae (30) 


so that on differentiating we find 
rd r(a — a 


If-in the last term of (30) we had used the effective area as intermediate 

between A —r and r we would have the factor (A — r) replaced by 

A —r/2. Thus to a first approximation we may ignore the factor 

(A — 2r/3)/(A — r) or (A — 2r/3)/(A — r/2) and again obtain 
r=7/Di . 

Thus m/(p, ¢) is given again by (24) for ¢ < #*. 


For ¢ > t* we have 


m* = my e-8Dt/ A*+3/4 | (32) 
so that 


m (p, t) = my e/te-suat (1-2). (337 


Now integrating m(p, t’) over t’ from zero to ¢ we obtain n(p, t) 


_ Mol PaCVy _ mo (A— p) Vit 
n (p, t) has tS aD? n (p, t) aN 
_ my (A— p)? (A-3p)? Ge ht mo (A— p) 34 
4dr) aD SS Zp (a) = 4) 
(7A +3 p) my (A— p) 3Dt/A2 A? 
SA hd) hel i! teak) kaa ary. AY 
3 SDs Seal 6) Moen Rares 


From expression (4), we may obtain the fraction of individuals who 
have been exposed to the information at any distance p from the center of 
the region. This fraction is given by 


N (p, &) aiae —mg(7A+3p)(A—p)/12D7N 


Nr Hh ; (35) 
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where NV, is the population density, taken to be independent of p: Let 
m = mo/ Nr, t* = t*/r. The total fraction of individuals exposed to the 
information is then 


2a ‘es : 1 AMA. 
ASAT pN(p wo) d = 1 —-—~ (1 — e7(7/8)m :* 
th? N 7 0 : e Mot * ( : 
ces (36) 
aie ae Lahn (5/3) V mot* 
—4vV mi* so ela) f eae rere 
(2/3) ¥ mot* 


If, for example, mol* K 1, » = nN A? is the total number of individuals 
in the region, if w is a width factor and is a probability of effective con- 
tact so that the mean free path \ = (pwN)-, expression (36) then be- 
comes 


—— iin Nrpwet .... (37) 


Hence the final value increases with m9 = mo/N7, the fraction started, 
with V,, the population density, and with the total population 7. 

Case 6. Continuous ‘‘point’’ sources of information: F < 1. One dimen- 
ston. Initially m = 0. Very rough approximation. Let S; be the strength of 
a point source of information, the source being in a large one-dimensional 
region. The information spreading out from the source will be dissipated 
as it travels since F < 1. The time integral of the source strength from 
zero to tire ¢ must be equal to the total amount which has moved away 
from the source due to diffusion and this in turn must equal the amount in 
the region at time ¢ plus that which has been lost due to F being less than 
one. As in the previous cases let the “‘particle’’ density be linear in the 
distance p from the source so that m(p, t) = m*(1 — p/r) where m* and 
r are functions of time. For p > r, m = 0. Then equating the derivatives 
of the above quantities, we have 


Dm dr dm* 1—F 
PL LS Te LL ae et 38 
hikes spss dae dt + pes ts as 
from which, if a = (1 — F)/r, 
aN nie (39) 
a 
Ne (40) 
a 


Therefore, after ¢ > 3/a,r ~ r* = VD/a = W711 — FB, m* = S/V/ Da. 
In order that r* be considerably larger than \, 1 — F must be very small. 
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Note that r* does not depend upon S considerably. In order to simplify 
matters we shall approximate (39) and (40) by the following: 


li vA 1 
p= De. m* = SVL, for tS) 


'D P S 1 
= 4\ se = —— ‘>=. 
r vee m rE for sar 


In this case if p = p/r*, i = t/r, N(p, #) is given by the following expres- 
sion, if¢ > 1/a or? > 1/(1 — F), 


N (p, t) = 1 — e~18(-)r1/18 1 F) 8/2) (8(1— F) 1-9} , (41) 


If t becomes large enough, it can be seen that 7(p, #) or V(p, #) approaches 
one if r < r*, whereas V(p, t) = 0 for r > r*. However, this latter result 
for r > r* is only an approximation so that although for some time there 
will be a moderately sharp boundary at r*, the value of WV will rise with 
time near r*, the rise being very slow for distances appreciably greater 
than 7*. But if a > 0, then the information is propagated as a wave 
spreading out over the region. 

Case 7. Informed individual may not respond to more than one hearing. 
Uniform spatial distribution. If an individual responds only to the first 
hearing, then expression (2) would be modified by replacing F by F times 
the fraction of individuals who have not yet been informed, 1 — V/N rp. 
Suppose now that a fraction # of the individuals do respond each time but 
the remaining fraction do not. Then instead of (2) we find 

om 


a 7 Vim + UF -1-F(1—p) NI™, (42) 


where V = N/Np. 
We shall consider the case in which m does not vary with position in 


order to relate the present method with that of Landau and Rapoport 
(loc. cit.). 


In this case we may write 
dm —_ m 


Introducing the value of W in terms of » as given by (4) and eliminating 
m by (3), we find, on integrating term by term, 


d 
1o= (Fp—1)n—F(1—p) Npe* T+ K, (44) 


where X is an integration constant. Since n/N above is the same — y of 
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equation (17) of Landau and Rapoport (Joc. cit.), it is evident that (44) 
above reduces to their expression (17) if = 0. Note that 7 and Fr are 
respectively 1/k and 4 in their expression (17). The only difference ap- 
pears in the meaning of the constant of integration. We have here sup- 
posed (0) and therefore V(0) = 0 so that the solution for p = 0 would be 
the same as that of equation (21) of Landau and Rapoport (1953), their z 
being equal to N, if zo in the lower limit of the integral is replaced by zero 


1,0 


Te 
No 


FIGuRE 1 


while z) under the integral is replaced by No, the value of NV (©) when 
F = 0. Since dn/dt = (dN/dt)/(1 — W), ¢ can be given as a function of 
N by a quadrature. If we wish the solution after a long time, we may set 
n = 0 in (44) and obtain an expression, which when rearranged may be 
written, if V¥. = W(@) (the fraction of those which are ultimately in- 
formed) 


(ae = ata ee 
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In Figure 1 are graphed values of V., as a function of No for p = 0, in 
which case expression (45) above is the same as expression (23) of Landau 
and Rapoport in which z* = WV, and z = No. In Figure 2 are graphed 
values of V., as a function of F for various values of p for the case in which 
im) — No > 0. If p = 1 the above expression (45) reduces to expression 
(10). Note that Vj is the same as 1 — e~” and is therefore the number of 


1,0 


fe) | 2 3 4 


-— 


FIGuRE 2 


individuals who are “started,” not counting individuals twice if they have 
accidentally received the information twice initially. ' 
The author is indebted to Drs. N. Rashevsky, A. Rapoport and H. G. 
Landau for reading the manuscript and for their many helpful sugges- 
tions. 
This work was aided in part by a grant from the Dr. Wallace C. and 
Clara A. Abbott Memorial Fund of the University of Chicago. 
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BOOK REVIEW 


Cyrit H. GouLpEN. Methods of Statistical Analysis. Second Edition, 1952. 462 pp. New York: 
John Wiley and Sons. $7.50. 


This book is written for agricultural research workers who already have some knowledge of 
statistical principles and want to learn the details of more advanced statistical calculations. 
Most of the book is taken up with worked-out examples giving the algebra and arithmetic of 
the statistical calculations involved in the design and interpretation of agricultural experi- 
ments. No knowledge of mathematics is assumed beyond elementary algebra, so that there is, 
of course, no attempt at derivations or proofs except for the simplest formulas. References to 
the literature are, however, given for the reader who wants to look further. 

The reason for this revision is to include new material incorporating the advances in statis- 
tics since the first edition of 1939. In this the author has been quite successful, bearing in mind 
the audience to which the work is addressed. The three chapters on design of (agricultural) 
experiments are especially thorough. There is a chapter on quality control and sampling inspec- 
tion, including sequential sampling, but the treatment is far too sketchy to be of use to the 
engineer. 

The book can be recommended for agricultural research workers; for others it may be of 
occasional use as a reference on the details of a statistical calculation. 


H. G. LAnpavu 
The University of Chicago 
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